Transcriptome Analysis Reveals the Resistance Mechanism of Pseudomonas aeruginosa to Tachyplesin I.

BACKGROUND: Tachyplesin I is a cationic antimicrobial peptide with a typical cyclic antiparallel β-sheet structure. We previously demonstrated that long-term continuous exposure to increased concentration of tachyplesin I can induce resistant Gram-negative bacteria. However, no significant information is available about the resistance mechanism of Pseudomonas aeruginosa (P. aeruginosa) to tachyplesin I.
MATERIALS AND METHODS: In this study, the global gene expression profiling of P. aeruginosa strain PA-99 and P. aeruginosa CGMCC1.2620 (PA1.2620) was conducted using transcriptome sequencing. For this purpose, outer membrane permeability and outer membrane proteins (OMPs) were further analyzed.
RESULTS: Transcriptome sequencing detected 672 upregulated and 787 downregulated genes, covering Clusters of Orthologous Groups (COGs) of P. aeruginosa strain PA-99 compared with PA1.2620. Totally, 749 differentially expressed genes (DEGs) were assigned to 98 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and among them, a two-component regulatory system, a beta-lactam resistance system, etc. were involved in some known genes resistant to drugs. Additionally, we further attempted to indicate whether the resistance mechanism of P. aeruginosa to tachyplesin I was associated with the changes of outer membrane permeability and OMPs.
CONCLUSION: Our results indicated that P. aeruginosa resistant to tachyplesin I was mainly related to reduced entry of tachyplesin I into the bacterial cell due to overexpression of efflux pump, in addition to a decrease of outer membrane permeability. Our findings were also validated by pathway enrichment analysis and quantitative reverse transcription polymerase chain reaction (RT-qPCR). This study may provide a promising guidance for understanding the resistance mechanism of P. aeruginosa to tachyplesin I.

Introduction

Pseudomonas aeruginosa (P. aeruginosa) is a Gram-negative and ubiquitous environmental bacterium. It is also an opportunistic human pathogen capable of causing a wide array of life-threatening acute and chronic infections, particularly in patients with compromised immune defense. A potential alternative to antibiotics for the control of P. aeruginosa and other microbes is antimicrobial peptides (AMPs). A number of AMPs have been reported to have a satisfactory inhibitory effect on P. aeruginosa, such as LL-37, tachyplesin, polymyxin, and helical peptides extracted from snakes.1,2

During co-evolution of hosts and bacterial pathogens, bacteria have developed an ability to sense and initiate an adaptive response to AMPs to resist their bactericidal activity. Previous studies indicated that Gram-negative pathogens exhibited resistance mechanisms of AMPs, including proteolytic degradation of AMPs, shielding of the bacterial surface, modification of the bacterial outer membrane, pumping AMPs into or out of the cell, downregulation of AMP expression,3,4 diverse genetic alterations, and alterations of amino acids and carbohydrate metabolism.5,6

Some bacteria (e.g., P. aeruginosa) also could develop resistance to a number of AMPs, such as polymyxin,7 and reports on resistance of polymyxins in P. aeruginosa have increased. A growing number of P. aeruginosa strains demonstrated resistance to AMPs due to mutations in two-component regulatory systems (e.g., PhoPQ and PmrAB). These two-component regulatory systems play substantial roles in lipid A modification systems, in which P. aeruginosa may become resistant to polymyxin.8 The arnBCADTEF operon encodes genes that modify lipopolysaccharide (LPS) by the addition of 4-aminoarabinose to the lipid A moiety, thereby reducing the negative charge on LPS and its ability to mediate the self-promoted uptake across the outer membrane of AMPs. The regulation of LPS modifications in P. aeruginosa is possibly involved in the complex regulatory network, mediating adaptive resistance to AMPs by at least 4 two-component regulatory systems (e.g., PmrAB, PhoPQ, ParRS, and CprRS).9 The molecular mechanisms of antibiotic resistance in P. aeruginosa are mainly involved in changes of porins and efflux pumps, however, it needs to be indicated whether the resistance of AMPs in P. aeruginosa is also involved in these changes?

Tachyplesin I, a 17-amino acid cationic antimicrobial peptide (CAMP), was originally isolated from hemocytes of marine horseshoe crabs in 1988.10 With potent and broad-spectrum activities against both Gram-positive and Gram-negative bacteria, tachyplesin I has been found as a promising candidate for the development of anti-infection, anti-tumor, and anti-virus drugs.11–13 Bacterial resistance to tachyplesin I has been induced in Aeromonas hydrophila XS91-4-1, P. aeruginosa CGMCC1.2620 (PA1.2620), P. aeruginosa ATCC27853 (PA27853), Escherichia coli ATCC25922 and F41 with cross-resistance to other AMPs (e.g., pexiganan, tachyplesin III, and polyphemusin I), conventional antimicrobial agents (amikacin), as well as the potential involvement of extracellular proteases in mediating tachyplesin I-resistance in Gram-negative bacteria reported in previous studies.14,15 A limited number of researches concentrated on the mechanism of bacterial resistance to tachyplesin I.14–16 However, the resistance mechanism of P. aeruginosa to tachyplesin I has largely remained elusive.

RNA sequencing (RNA-Seq) is an advanced method for transcriptome profiling that can provide more precise measurement for genome-wide quantitative analysis of all types of transcripts (including coding, non-coding, and antisense) with a relatively high accuracy. RNA-Seq can reveal specific biological processes, such as biological evolution, basic metabolism, disease occurrence, and drug targets. In addition to the extensive utilization of RNA-Seq in clinical researches, it is very sensitive for measuring low abundance genes.17 In the present study, we employed RNA-Seq to investigate differences in genome-wide gene expression between tachyplesin I-resistant PA-99 and the original strain PA1.2620, followed by validation of our findings via pathway enrichment analysis and quantitative reverse transcription polymerase chain reaction (RT-qPCR). Then, the bacterial outer membrane permeability and outer membrane proteins (OMPs) were further analyzed to reveal the resistance mechanism of P. aeruginosa to tachyplesin I. The findings may contribute to enhance our understanding of the resistance mechanism of P. aeruginosa to tachyplesin I.

Materials and Methods

Strains, Media, and Growth Conditions

PA1.2620 was provided by China General Microbiological Culture Collection Center (CGMCC; Beijing, China). P. aeruginosa ATCC27853 (PA27853) was provided by the laboratory of the College of Life Science and Engineering, Henan University of Urban Construction (Pingdingshan, China). The minimal inhibitory concentration (MIC) values of tachyplesin I were 10 and 5 μg/mL for PA1.2620 and PA27853, respectively, in which these two strains were sensitive to tachyplesin I. Tachyplesin I-resistant mutants (PA-99 and P. aeruginosa ATCC27853-88-2 (PA-88)) were synthesized in our laboratory as described previously.14,15 These strains were cultured in Mueller-Hinton broth (MHB) medium, and then seeded onto nutrient agar plates at 30°C or 37°C.

Tachyplesin I

Tachyplesin I (>95% purity) was synthesized by Gil Biochemical Co., Ltd. (Shanghai, China), and its sequence was as follows: NH2-K-W-C-F-R-V-C-Y-R-G-I-C-Y-R-R-C-R-CONH2, including two disulfide bonds (C3–C16 and C7–C12).10 It was dissolved in sterile water to yield a 10 mg/mL stock solution, which was filter-sterilized before utilization. Peptide solution was freshly prepared on the day of the assay or stored at −20°C for a short period of time.

Selection of PA-99 and PA-88 Mutant Strain Resistance to Tachyplesin I

Bacteria were cultured in MHB with constant shaking at 160 rpm at 30°C (PA1.2620) or 37°C (PA27853). We transferred bacteria daily by inoculating 20 μL stationary phase culture in 2 mL MHB. All cells were initially grown in medium without tachyplesin I for five transfers. At the 6th transfer, 20 μL cell suspension was added to 2 mL MHB with or without tachyplesin I at a final concentration equal to half of MIC with shaking for 15 transfers at 160 rpm for 20 h. The bacteria were thereafter transferred into broth containing a double concentration of tachyplesin I every 10 transfers or more frequently, if the selected strains showed weak growth. The experiment was conducted for 88–99 serial transfers. All experiments were performed in duplicate.

RNA Isolation, cDNA Library Construction, and Transcriptome Sequencing

PA1.2620 original strain and PA-99 mutant were grown overnight for within 12 h in 3 mL MHB medium at 30°C, and shaken at 180 rpm. Then, cultures were refreshed with medium to OD600nm=0.2 and grown to mid-exponential phase (OD600nm=1.0) until the time of harvesting. Extraction and purification of the total RNA from the cells were carried out using a Bacterial Total RNA Purification Kit (Sangon Biotech Co., Ltd., Shanghai, China). To ensure the use of total RNA sequencing for transcriptome analysis, quality, quantity, and integrity of extraction were analyzed by using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), gel electrophoresis, and an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA), respectively. All experiments were performed in triplicate.

After the samples of total RNA passed the test, the library construction was carried out. The main procedures for library construction were as follows: First, remove 16S ribosomal RNA (16S rRNA) of prokaryotes with a probe complementary to the rRNA gene sequencing technique. Second, join the fragmentation buffer in rRNA-depleted RNA and randomly interrupt it. Third, deplete rRNA as a template, synthesize the first cDNA chain using six random primers, then, add the buffer, dUTP, dATP, dGTP, dCUP, RNase H, and DNA polymerase I to synthesize the second cDNA chain using AMPure XP Beads purifying cDNA. AMPure XP Beads were utilized to select strand size. Special enzymes were used to degrade the second strand of cDNA containing U bases. Finally, the cDNA library was constructed by PCR enrichment.

The RNA-Seq libraries of six samples were prepared and subsequently sequenced using the Illumina HiSeq 2500 Platform (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s instructions. The construction of the libraries and the transcriptome sequencing were performed by Biomarker Biotechnology Co., Ltd. (Beijing, China).

Analysis and Assembly of RNA-Seq Data

Based on sequencing by synthesis technology, the clean reads were assembled using the Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA). Generally, the raw reads contain tiny minority primer sequences, adapter sequences, sequencing connector, and other potential contaminants. Prior to subsequent analysis, the clean reads were filtered from the raw reads by removing the reads with only adaptor and unknown nucleotides.

The clean reads obtained from RNA-Seq were mapped on the reference genome of P. aeruginosa PAO1 by Bowtie-2 software, and only mapped data were used for subsequent analysis.18 Data analysis and base-calling were undertaken by using Illumina sequencing software (Illumina Inc., San Diego, CA, USA).

Gene Expression Analysis

To evaluate the expression level of all reads, each UniGene cluster was compared via UniGene database by using Bowtie-2 software, and normalized into RPKM (Reads Per Kilobase of transcript per Million fragments mapped) values according to the following formula: RPKM = (Map the UniGene cluster)/[Map all UniGene clusters (million)]) × (The length of UniGene cluster (kb)).1 Besides, differentially expressed genes (DEGs) with a false discovery rate (FDR < 0.01) and a fold change of ≥2 were selected for analysis.19 Functional annotation of each UniGene cluster was searched against various protein databases, and identified by annotation information of the given UniGene cluster that has the highest sequence similarity with the tested one. All UniGene clusters were also searched in the National Center for Biotechnology Information (NCBI)’s non-redundant protein database (Nr), Swiss-Prot database, and Cluster of Orthologous Groups (COGs) of proteins using BLAXTX with an E-value of less than 10–5.20 The Blast2GO program and WEGO software were also employed to obtain the gene ontology (GO) annotations and the distribution of gene functions for each UniGene cluster using a value of less than 10–5.21 Assignment of each UniGene cluster to different pathways was carried out by searching in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using the KEGG Automatic Annotation Server.22 DEGs were annotated in the Comprehensive Antibiotic Resistance Database (CARD) to find out the antibiotic resistance genes.23

Validation by RT-qPCR

To increase the reliability of the analysis, RT-qPCR was conducted to verify some of the changes in expressions detected by RNA-Seq. Herein, 10 target genes were respectively selected from PA-99 mutant and the original strain PA1.2620 for analysis by RT-qPCR. For two strains, RT-qPCR was carried out in triplicate. Moreover, RT-qPCR was undertaken on an ABI 7900 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) using KAPA SYBR FAST qPCR kit (Kapa Biosystems Inc., Wilmington, MA, USA). Specific primers of the genes showed in were designed using Primer Premier 5 software based on the genome sequence of P. aeruginosa. The thermal cycler conditions were as follows: 40 cycles of denaturation at 94°C for 20 s, annealing at 60°C for 20 s, and extension at 72°C for 30 s. Using 16SrRNA as a reference, relative expression levels of the genes were quantified using 2−ΔΔCt method.24

NPN Assay to Determine Outer Membrane Permeability

Outer membrane permeability was determined by using the NPN assay with minor modification.25 Cells were harvested from overnight cultures by centrifugation, washing, and re-suspending in 25 mM Tris-HCl (pH 7.2). The value of OD600 was adjusted to 0.5 with the same buffer; however, the exact OD600 was measured and recorded. In addition, 80 μL of NPN (1 mM) was swiftly mixed with 1.92 mL of the above-mentioned cell suspension. The fluorescence of the mixture was monitored immediately by a spectrofluorometer (PerkinElmer, Waltham, MA, USA). The wavelengths of excitation, emission, and slits were 350, 420, and 5 nm, respectively. The permeability was indicated by the fluorescence absorption per OD600 value of the sample.

Analysis of OMPs

Preparation of OMPs and Determination of Their Total Contents

OMPs of the original and mutant strains were extracted by a previously described method26 with minor modification. A single colony was selected, which inoculated in 2 mL MHB medium, P. aeruginosa strains at 30°C or 37°C for 18 h at 180 r/min. After the culture was inoculated in 100 mL MHB medium for 4 h under the same condition, the cells were centrifuged (at 6000 ×g for 10 min at 4°C), suspended in 15 mL of 10 mM phosphate-buffered saline (PBS, pH 7.0) containing 0.14 mol/L β-mercapto ethanol, and disrupted by ultrasonic wave for 6 min in ice bath. Cellular debris was removed by centrifugation at 6000 ×g for 10 min at 4°C. The supernatant was centrifuged at 10,000 ×g for 30 min at 4°C. Outer membrane precipitation collected was added to 5 mL of 2% sodium lauroyl sarcosinate (Sarkosyl; Sigma-Aldrich, St. Louis, MO, USA) and 5 mL of 10 mM PBS (pH 7.0) incubated at room temperature overnight. The mixture was centrifuged at 20,000 ×g for 90 min to recover the detergent-solubilized OMPs. The pelleted proteins were resuspended in 3 mL of 10 mM PBS (pH 7.0), incubated with 1% sarkosyl at room temperature for 20 min, and recovered by ultracentrifugation as described above. The final pellet was re-suspended in 1 mL of 10 mM PBS and stored at −20°C until analysis.

The content of OMPs was determined by ultraviolet absorption spectroscopy according to the following formula: protein concentration (mg/mL) = 1.45 * OD280 - 0.74 * OD260. Protein concentration was adjusted for a final concentration of 45 μg/mL, and the protein samples were stored at –80°C.

SDS-PAGE Assay and Calculating the Molecular Weight of the Protein

Sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) was carried out with 12% separating gel and 5% concentrating gel after preparation of OMPs was solubilized according to 1:4 (V/V) at 100°C for 5 min in 1 M Tris-HC1 buffer (10% SDS, 5% 2-mercaptoethanol, 50% glycerol, and 0.1% bromophenol blue). When the molecular weight of the protein was between 15,000 and 200,000 Da, electrophoretic mobility and logarithm of molecular weight showed a linear relationship according to the equation of LgMr = K - bmR; where Mr represents the molecular weight of protein; K is a constant; and mR denotes the relative mobility. Under certain conditions, both b and K are constant. A standard curve can be obtained by plotting the logarithm of the molecular weight of standard proteins with known molecular weight. Those proteins with unknown molecular weight underwent electrophoresis under the same conditions, in which the molecular weight could be obtained by the standard curve according to its electrophoretic mobility.

Results

Selection of Tachyplesin I-Resistant PA-99 and PA-88 Mutants

We experimentally induced bacterial resistance to tachyplesin I in P. aeruginosa in vitro for long term, in which concentration of tachyplesin I continuously increased. High-resistance to tachyplesin I was developed after 99 serial transfers in PA1.2620 and 88 serial transfers in PA27853. There were significant differences in the MIC values for PA1.2620 and PA27853 between the induced strains and the corresponding induced strains for treatment. The MIC values of tachyplesin I were 150 and 80 μg/mL for PA-99 and PA-88 mutants, respectively, in which these two mutants were highly resistant to tachyplesin I. Considering the expenses of the experiments, the similarity of strain, and the higher resistant to tachyplesin I, we only investigated differences in genome-wide gene expression between tachyplesin I-resistant PA-99 and PA1.2620 by RNA-Seq analysis.

Sequence Analysis and Assembly

Transcriptome sequencing of six cDNA libraries was performed by using total RNA extracted from tachyplesin I-resistant PA-99 strain and the original strain PA 1.2620. After filtering raw data and passing quality control, clean data of each sample were no less than 2.11 Gb, GC content of clean data was above 62.57%, and Q30 base reached over 95.42% (Table 1). These data were selected as high-quality reads for further analysis.

Table 1

Transcriptome Sequencing Summary for Different P. aeruginosa Strains

#Sample ID	ReadSum	BaseSum	GC (%)	Q30 (%)
PA1.26201	12,794,030	3,699,522,416	64.49	96.43
PA1.26202	14,442,212	4,311,426,684	63.34	95.46
PA1.26203	12,387,552	3,694,243,854	63.73	96.21
PA-991	10,524,519	3,139,186,296	63.65	96.32
PA-992	10,602,206	3,161,091,324	62.57	95.60
PA-993	11,337,636	3,384,373,320	63.54	95.42

Notes: Superscripted 1-3 stands for three biological duplications of one sample. BMK_ID: the uniform number of the sample from BMCG; ReadSum: total number of pair-end Reads in Clean Data; BaseSum: clean data total base number; GC (%): Clean data GC content, that is, the percentage of G and C bases in Clean Data in total bases; Q30 (%): the percentage of bases with a clean data mass value greater than or equal to 30.

By evaluating the transcriptome of sequencing data of six samples, including those mapped to the reference genome, randomness testing of mRNA fragmentation, saturation testing of data, and distribution of reads on the reference genome, the results revealed that the transcriptome data obtained were ideal and reliable in terms of quality. The data related to quality and mapping statistical data are presented in .

In the result of transcriptome sequencing, we found a correlation between one biological replicate and the other two replicates, with low correlations (0.4689, 0.4735) among the three biological replicates of the original strain PA1.2620, whereas R2 between the other two samples was 0.9931 (see ). The reason may be that the one biological replicate was prepared and tested without any association with the other two biological replicates. Considering the accuracy of the results, we deleted one biological replicate with a low correlation. Therefore, the expression profiling analysis was undertaken based on the two biological replicates of PA1.2620.

Gene Expression Analysis and Functional Annotation

To investigate the changes in transcription levels before and after PA1.2620 resistant to tachyplesin I, we selected tachyplesin I-resistant PA-99 mutant and the original strain PA1.2620 to perform expression profiling analysis. Results of hierarchical clustering and volcano plots are shown in Figure 1A and B. As illustrated in Figure 1A, PA-99 resistant mutant had several DEGs compared with the original strain, and log2FC values of expression levels were mainly distributed at the range of −5 and 5. There were 1459 DEGs with 318 novel genes, including 672 (10.64%) upregulated genes and 787 (12.45%) downregulated genes among 6318 genes. There were 254 downregulated genes in the 318 novel genes. Besides, the expression of 57 genes was upregulated with more than 5-fold changes, and the expression of 30 genes was downregulated with more than 20-fold changes. According to BLAST analysis, the three upregulated and downregulated genes were described as hypothetical proteins in the DEGs. The novel gene (gene ID: Novel_392) and PA2754a were both downregulated form of the above-mentioned genes, and no homologous proteins were found by BLAST analysis. The third downregulated gene was PA2764 (log2FC=−16.39), and its homologous protein was alpha/beta hydrolase superfamily. The most significantly upregulated gene was PA0912, which had high homology with the lipoprotein-related gene of Pseudomonas stutzeri. The second upregulated gene was Novel_146. The third upregulated gene was Novel_94 (Novel_94), and its homologous protein was RNA polymerase-binding transcription factor (DksA). It needs to be further indicated whether PA-99 resistant to tachyplesin I was associated with the changes in highly expressed genes described as hypothetical proteins. Numerous DEGs related to drug resistance were identified, which mainly involved in efflux pumps, porins, OMPs, two-component regulatory systems, etc. We speculated that these DEGs may play a substantial role in PA-99 resistant to tachyplesin I.

Figure 1

Hierarchical clustering and Volcano plot of differential expression in PA1.2620 original strain vs PA-99 mutant. (A). Hierarchical clustering. T04-T05 stands for two biological duplications of PA1.2620 strain, T09-T11 stands for three biological duplications of tachyplesin I-resistant PA-99 strain. Each column represents one sample, and each row represents one gene. The red and green gradients indicate an increase and a decrease in gene expression abundance, respectively. (B). Volcano plot. Each point in the volcano plot represents one gene, and the abscissa represents the logarithm of the expression difference fold of one gene in PA1.2620 original strain vs PA-99 mutant. This vertical axis represents the negative logarithm of the error detection rate. The green dots represent downregulated genes, the red dots represent upregulated genes and the black dots represent unchanged genes.

All DEGs of PA1.2620 versus PA-99 were tested for functional annotation using COG, GO, KEGG, Swiss-Prot, NR, and Pfam databases, in which the number of DEGs annotated in different databases is presented in Table 2.

Table 2

The Numbers Statistics of Annotated DEGs

Type	COG	GO	KEGG	Swiss-Prot	NR	Pfam
PA1.2620 vs PA-99	960	1,004	749	979	1,424	1,198

A further functional classification of all DEGs in 5 samples was performed using the GO database. The DEGs of GO annotation classification indicated that P. aeruginosa resistant to tachyplesin I was enriched in cellular components (CCs), molecular functions (MFs), and biological processes (BPs), which was subdivided into 43 functional groups, as shown in Figure 2A and . In BPs, DEGs were mainly related to localization, single-organism process, cellular process, and cell killing. In CCs, the differences in enrichment trends were related to cell (26.49% vs 31.27%), macromolecular complex (3.78% vs 4.75%), organelle (1.59% vs 2.67%), and cell part (26.49% vs 31.27%) in DEGs. For MFs, the differences in enrichment trends in DEGs were related to receptor activity (1.49% vs 2.44%), transporter activity (11.35% vs 12.69%), and binding (46.51% vs 47.64%). Collected data from CCs, MFs, and BPs suggested that P. aeruginosa resistant to tachyplesin significantly affected cell, cell part, receptor activity, transporter activity, and binding.

Figure 2

Functional annotation of DEGs in PA1.2620 original strain vs tachyplesin-resistant PA-99 mutant. (A). GO enrichment categorization; (B). COG classification; (C). KEEG categorization.

Next, DEGs were further analyzed by COG enrichment analysis. COG was classified into 21 functional categories, as illustrated in Figure 2B and . The top 3 COG classes based on the percentage in total COG were as follows: general functional prediction (20.20%), amino acid transport and metabolism (15.10%), and transcription (13.54%), in which they totally accounted for 48.85% of COG. These 3 COG classes may be closely associated with P. aeruginosa resistant to tachyplesin I. The following ranked COG was unknown function (10.52%), inorganic ion transport and metabolism (9.48%), signal transduction pathways (9.28%), energy production and conversion (8.85%), cell wall/membrane/envelope biogenesis (7.19%), posttranslational modification, protein turnover, chaperones (6.67%), carbohydrate transport and metabolism (6.67%), and lipid transport and metabolism (6.56%). In the above-mentioned 11 COG categories, each category accounted for over 6%, which may be related to resistance to tachyplesin I. Meanwhile, COG categories without functional changes included RNA processing and modification, nuclear structure, as well as cytoskeleton and extracellular structures, indicating that P. aeruginosa resistant to tachyplesin I did not have any correlation with these functions.

In total, with the help of KEGG pathway analysis, 749 DEGs were assigned to 98 KEGG pathways. These pathways were mainly divided into “the genetic information processing”, “cellular processes”, “environmental information processing”, “human diseases”, and “metabolism” (Figure 2C). Among them, the metabolism pathway was involved in the majority of DEGs, including 40 DEGs with the biosynthesis of amino acids (17 upregulated and 23 downregulated genes), carbon metabolism (41 DEGs), pyruvate metabolism (22 DEGs), purine metabolism (20 DEGs), etc. These results showed that the tachyplesin I-resistant PA-99 mutant resulted in the alteration of several metabolic pathways. The two-component regulatory system (ko02020), beta-lactam resistance system (ko01501), lipopolysaccharide biosynthesis (ko00540), and ABC transporter (ko02010) pathways were involved in a number of known drugs resistant to DEGs. In addition, 43 DEGs with 18 novel genes included 23 upregulated and 20 downregulated genes in a two-component regulatory system, in which these genes were mainly related to OmpR family (PhoB, PhoA, PhoP, OprD, etc.), NarL family, NtrC family, and Chemotaxis family (). There were 36 DEGs, including 21 upregulated and 15 downregulated genes in ABC transporter (). Additionally, 6 DEGs (3 upregulated and 3 downregulated genes) in beta-lactam resistance system were involved in genes encoded by outer membrane porin oprD (gene ID: Novel_156 and Novel_157), efflux pumps of the resistance-nodulation-division family (mexB, PA4374), mexR for repressor of the mexAB-oprM, and armR anti-repressor for mexR ().

Annotation and Analysis of Resistant Genes

Firstly, to select the known resistant genes, DEGs were annotated in the CARD database, and a total of 14 antibiotic resistance genes were matched. The resistance mechanisms of these genes mainly involved in antibiotic target alteration (pmrB, mexR), antibiotic efflux (PA2018, PA0156, PA4136, PA3677, PA4374, mexC, mexE, mexB, armR, nalC, PA3521) mainly coding RND efflux pumps, and antibiotic inactivation (PA5514 coding for OXA beta-lactamase). The majority of these genes were upregulated, indicating that the resistance mechanism of P. aeruginosa to tachyplesin I may be related to overexpression of efflux pumps.

Secondly, except for the selected antibiotic resistance genes, according to the gene annotation, DEGs contributed to resistance to antibiotic or AMPs were classified and statistically analyzed. These genes were mainly related to encoded efflux pumps, a two-component regulatory system, a DNA mismatch repair system, etc. (Table 3). Besides, 8 genes from the major facilitator superfamily (MFS) were detected as antibiotic resistance genes, such as PA1131, PA0334, and PA1282 encoded by tetracycline resistance protein, PA4136 encoded by bicyclomycin resistance protein, and PA3573 encoded by multidrug resistance transporter. Moreover, 4 genes (PA2837, PA3404, PA3521/opmE, and PA4592) encoded by outer membrane efflux protein were all upregulated (log2FC≥1.236), as well as 7 genes encoded by outer membrane porin, oprD family with 3 novel-genes and 4 known genes; except for opdB that was upregulated, others were all downregulated. In particular, the PA2760 (oprQ) gene was strongly downregulated (log2FC≥14).

Table 3

Expression of Part Genes Related to Drug Resistance

	Gene Name or #ID	Gene Description	PA1.2620 vs PA-99
FDR	log2FC	Regulated
RND efflux pump and regulator gene	PA0156 (triA)	RND efflux membrane fusion protein	0.0092193	1.232266998	Up
	PA1435	RND efflux membrane fusion protein	1.82E-05	1.241591925	Up
	PA1436	RND efflux transporter	–	–
	PA3676 (mexK)	RND efflux transporter	–	–
	PA3677 (mexJ)	RND efflux membrane fusion protein	0.000985076	1.40757487	Up
	mexA	Multidrug resistance membrane fusion protein MexA	–	–	–
	mexB	Multidrug resistance protein MexB	2.50E-11	1.190364744	Up
	oprM	Outer membrane protein OprM	–	–
	PA4374 (mexV)	RND efflux membrane fusion protein	0.00010801	1.612935679	Up
	mexC	RND multidrug efflux membrane fusion protein MexC	1.95E-10	1.936213338	Up
	mexD	RND multidrug efflux transporter MexD	–	–
	oprJ	Multidrug efflux outer membrane protein OprJ	–	–
	mexE	RND multidrug efflux membrane fusion protein MexE	4.86E-06	1.314524548	Up
	mexF	RND multidrug efflux transporter MexF	–	–
	oprN	Multidrug efflux outer membrane protein OprN	–	–
	mexH	RND efflux membrane fusion protein	–	–
	mexI	RND efflux transporter	–	–
	opmD	Hypothetical protein	–	–
	PA3521 (opmE)	Outer membrane efflux protein	0.000128	1.545374	Up
	PA2837	Outer membrane efflux protein	3.35E-10	1.458729	Up
	mexR	Repressor of the mexAB-oprM multidrug resistance operon	0.000572908	−1.286050369	Down
	armR	Antirepressor for MexR	0.0034723	1.311197	Up
nalC	RND transcriptional regulator	5.05E-09	−1.0027	Down
Two-component regulator system	phoB	Two-component response regulator PhoB	0.000531045	1.178656488	Up
	pmrB	Two-component regulator system signal sensor kinase pmrB	2.27E-05	1.000160392	Up
	pmrA	Two-component regulator system response regulator PmrA	1.23E-18	1.649392693	Up
	Novel_194	Response regulator PhoP	0.001111	1.924075	Up
	Novel_498	Alkaline phosphatase (PhoA)	0.000897	2.476633	Up
	Novel_499	Alkaline phosphatase (PhoA)	0.000995	2.513525	Up
	PA0034	Multidrug resistance, response regulator EvgA	4.68E-06	−1.14216	Down
Novel_346	Sensor histidine kinase EvgS	0.000589	−1.13806	Down
MFS resistance protein	PA3137	Multidrug export protein EmrB	0.001280967	−1.186834062	Down
	bacA	Bacitracin resistance protein BacA	2.13E-32	−3.29214	Down
	PA1282	Antiseptic resistance protein	1.74E-21	−2.13973	Down
	PA1316	Multidrug resistance protein	4.36E-13	−1.70267	Down
	PA4136	Bicyclomycin resistance protein	0.00207	−1.316	Down
	PA0334	Tetracycline resistance protein signature	1.30E-11	−1.5386	Down
	PA1131	Tetracycline resistance protein	7.88E-05	1.05043	Up
PA3573	Multidrug resistance transporter, Bcr/CflA family	2.56E-05	2.067501304	Up
Outer membrane efflux protein	PA3404	Outer membrane efflux protein	8.72E-06	1.23575	Up
PA4592	Outer membrane efflux protein	0.000174	1.496072	Up
Outer membrane porin D, OprD family	opdB	Outer membrane porin D,	0.009088	1.05881	Up
	Novel_156	Outer membrane porin D	1.46E-11	−2.66032	Down
	Novel_157	outer membrane porin D	1.17E-06	−2.014076	Down
	PA2760 (oprQ)	outer membrane porin D	1.46E-155	−14.4306	Down
	PA4179	outer membrane porin, OprD family	0.001555	−1.1668	Down
	Novel_314	Outer membrane porin D	0.000316	−2.26031	Down
Novel_443	Outer membrane porin D	9.40E-06	−2.09792	Down
Outer membrane protein	opr86	Outer membrane protein Opr86	0.000113695	−1.514177711	Down
lptF	Outer membrane protein OmpA	2.47E-07	−1.24917	Down
DNA mismatch repair protein	mutL	DNA mismatch repair protein	8.83E-07	1.085268	Up
	mutS	DNA mismatch repair protein MutS	0.005549	1.130611	Up
	radA	DNA repair protein RadA	2.62E-06	−1.32025	Down
radC	DNA repair protein RadC	6.65E-06	−1.42294	Down
Cell wall/membrane/envelope biogenesis	arnT	Cell wall/membrane/envelope biogenesis	0.001764	1.162594	Up
	murD	UDP-N-acetylmuramoylalanine-D-glutamate ligase	7.06E-05	1.475058	Up
	waaC	Lipopolysaccharide heptosyltransferase	8.06E-05	1.023263	Up
lasA	Protease LasA	2.62E-07	−1.73674	Down
beta-lactamase	PA2315	Beta-lactamase	5.16E-05	1.92727	Up
PA5514	Beta-lactamase	1.41E-06	−1.41767	down
Biofilm formation	chpC	Chemosensory pili system protein ChpC/chemotaxis protein	0.004917388	1.704317455	Up

According to the above-mentioned results, RND, MFS, and oprD families, and a two-component regulatory system play a pivotal role in P. aeruginosa resistant to tachyplesin I.

Validation via RT-qPCR

To validate the results of transcriptome sequencing, 10 DEGs from tachyplesin I-resistant mutant PA-99 and original strains (e.g., PA1435, PA2018, mexB, PA4396, parR, phoB, PA3573, PA0156, mexC and motB) were randomly selected for RT-qPCR. As displayed in Figure 3, all the selected DEGs showed concordant expression patterns between the RNA-Seq and the results of RT-qPCR.

Figure 3

Illustrating the RT-qPCR conformation results for the 10 selected DEGs. The X-axis represents the selected 10 genes and the Y-axis represents the fold change values derived from RNA-seq and RT-qPCR.

NPN Assays

To indicate whether P. aeruginosa resistant to tachyplesin I was involved in the changes of the outer membrane permeability, an NPN (n-phenyl-1-naphthylamine) assay was performed in different P. aeruginosa strains. The nonpolar probe NPN tends to stay in a hydrophobic environment and fluoresces strongly in a hydrophobic environment, while weakly in an aqueous environment.

The results of this assay revealed that two tested mutants, namely PA-99 and PA-88 mutants, exhibited weaker outer membrane permeability compared with PA 1.2620 and PA27853, respectively (Figure 4). This provides a tenable explanation for resistance to tachyplesin I, as the peptide would be able to translocate across the outer membrane more difficult than in the original strain. These data indicated that P. aeruginosa resistant to tachyplesin I was associated with the decrease of outer membrane permeability.

Figure 4

Outer membrane permeability of original and mutant strains.

To clarify whether OMPs of P. aeruginosa play a role in outer membrane permeability and drug resistance, we detected the content and differences in OMPs in two strains by SDS-PAGE. As shown in Table 4, the total concentration of OMPs was higher in both original strains than that of the resistant-tachyplesin mutants. Figure 5 illustrates the molecular weight and the amount (%) of OMPs extracted from original strains using SDS-PAGE. Additionally, the bands of OMPs were similar in both original strains (Figure 5). According to the standard curve, the relative molecular weight of the tested protein was obtained. A linear regression equation was achieved according to the standard curve: y=−1.0858x+5.4353, R2=0.9934. As depicted in Figure 5A, there were three bands with molecular weights of 69.28, 43, and 46.19 kDa between PA1.2620 original strain and tachyplesin I-resistant PA-99 strain. The molecular weights of two bands (69.28 (A) and 43 kDa) were reduced, while the molecular weight of one band (46.19 kDa (C)) was increased in PA-99 mutant. Similarly, compared with PA27853, the molecular weight of one band (69.28 (A) kDa) was attenuated, while molecular weights of two bands (46.19 (C) and 52.87 kDa (B)) were increased in PA-88 mutant (Figure 5B). The results showed that OMPs in both tachyplesin-resistant P. aeruginosa strains were different.

Table 4

The Contents of OMPs from Different P. aeruginosa Strains

Strains	Protein Concentration (mg/mL)
PA1.2620	0.92±0.25
PA-99	0.66±0.21
PA27853	0.69±0.15
PA-88	0.48±0.18

Figure 5

Analysis of OMPs of different P. aeruginosa strains. (A). 1, 4: PA1.2620; 2–3: PA-99. (B). 1, 4: PA27853; 2–3: PA-88.

Abbreviation: M, standard protein marker.

Discussion

To our knowledge, this research is the first attempt to explore the resistance mechanism of P. aeruginosa to tachyplesin I using whole-genome expression profiling. For detecting the changes of gene expression in tachyplesin-resistant PA-99 strain relative to the original PA1.2620 strain, transcriptome sequencing was carried out, and the results were further validated by RT-qPCR. From the transcriptome data, we found that a great number of genes exhibited significantly altered expressions.

In order to enhance the accuracy and reliability of the analysis, each test was performed in triplicate in our study. With respect to Pearson’s correlation coefficient (R) as an indicator for evaluation, R2 is the closer to 1, and the correlation between the two repeated samples is the stronger.27 Therefore, in the present study, considering the accuracy of the results, we deleted one biological replicate in the original PA1.2620 strain, which was poorly correlated with the other two samples ().

A two-component regulatory system plays a substantial role in the pathogenicity, virulence, biofilm formation, and drug resistance in P. aeruginosa. In the current research, we detected 43 DEGs in a two-component regulatory system, including some regulated resistant genes, such as pmrA-pmrB, response regulator phoP (gene ID: Novel_194, log2FC=1.924075), alkaline phosphatase encoded by phoA (gene ID: Novel_498, 499), and a two-component response regulator phoB in PA-99 mutant. A two-component regulatory system may also play a significant role in P. aeruginosa resistant to tachyplesin in the present study. A previous research demonstrated the importance of regulatory genes, such as phoP-phoQ and pmrA-pmrB, as well as the addition of aminoarabinose to Lipid A.8 The pmrB mutant of P. aeruginosa showed a significant metabolic perturbation in the total intracellular lipid level, the methionine salvage cycle, and synthesis of spermidine.5 Resistance can be mediated by the addition of positively charged arabinosamine through the action of the arnBCADTEF operon. The results achieved in the present research were consistent with a previous report, in which the complex regulation of LPS modification involved the participation of at least 4 two-component regulatory systems in P. aeruginosa.9

Drug efflux through efflux pumps is one of the main resistance mechanisms of antibiotic or AMPs of P. aeruginosa. Some bacteria can use efflux pumps or increase the expression of efflux pumps to mediate resistance against CAMPs.3,28 Especially, the efflux pumps of the RND family were studied in P. aeruginosa, which could be directly related to the efflux of antibiotics. At least 7 RND family drug efflux pumps have been reported to exist in cells of P. aeruginosa, such as MexAB-OprM, MexCD-OprJ, MexEF-OprN, MexXY-OprM, MexJK-OprM, MexVW-OprM, and MexGHI-OpmD, which could be involved in the antibiotic resistance.29 MFS is another important superfamily of secondary membrane transporter, playing a significant role in substance exchange, energy metabolism, and drug resistance. The ABC family includes transporters involved in both uptake and efflux, and these transporters can transport a wide range of substrates, including sugars, amino acids, ions, drugs, polysaccharides, and proteins. In the present study, DEGs were found in ABC transporters, MFS transporters, and RND superfamily. There were a number of upregulated genes, which were encoded by RND efflux membrane fusion proteins and outer membrane efflux proteins. For instance, there were 7 upregulated genes in RND superfamily, including 6 efflux membrane fusion proteins (mexC/E/V/J, PA0156, PA1435), 1 export membrane protein (mexB), and 1 downregulated gene that is responsible for repressor of the mexAB-oprM multidrug resistance operon (mexR). No encoded multidrug efflux outer membrane proteins (oprJ/M/N, opmD) were expressed, while genes with unknown function (opmE, PA2837), which were encoded by outer membrane efflux proteins, were upregulated (Table 3). Besides, 31 DEGs of MFS, including 15 up- and 16 downregulated genes, contained 8 known genes related to antibiotic resistance. In ABC transporters, there were 36 DEGs; however, whether they are related to the resistance mechanism of P. aeruginosa should be further clarified. The expression levels of efflux pumps showed that P. aeruginosa resistant to tachyplesin was also related to the expression of efflux pumps, while their participation in drug resistance needs to be further studied.

Another mechanism of bacterial resistance to AMPs is the reduced entry of the AMPs into the bacterial cell due to modifications of the cell surface that limit interactions with the drug or reduce the number of entry channels (porins). OMPs of P. aeruginosa play a key role in membrane permeability, antibiotic resistance, nutrient uptake, and virulence in the infection site. In addition, reduction or loss of the oprD porin can reduce or block the antibacterial drug to enter into P. aeruginosa in part due to the decreased permeability of the P. aeruginosa outer membrane.30 To investigate whether P. aeruginosa resistant to tachyplesin I could be involved in expression of the genes encoded by OMPs or changes of outer membrane permeability, transcriptional regulation revealed that 4 genes encoded by outer membrane efflux protein were all upregulated (log2FC≥1.236 folds), 6 genes encoded by outer membrane porin D were downregulated, and only 1 opdB gene was upregulated in resistant strains of PA 99 (Table 3). The genes (ID: Novel_156 and Novel_157) encoded by oprD belonged to β-Lactam resistance and two-component regulatory systems. In particular, the PA2760 (oprQ) gene encoded by porin was strongly downregulated (log2FC≥14). The previous study indicated that P. aeruginosa oprQ protein did not function in roles common to this class of OMPs. The regulation of oprQ expression may be in dual regulation with other proteins and the overexpression of oprQ could lead to increased sensitivity to certain classes of antibiotics.31 A co-regulation of the two porin genes (oprD and oprQ) has been recently reported since both genes are overexpressed in an Hfq mutant of P. aeruginosa.32 We further analyzed the changes of outer membrane permeability and OMPs, and the results showed that the reduced outer membrane permeability and decreased levels of OMPs were related to tachyplesin I-resistant P. aeruginosa strain compared with the original strain. According to the achieved results, it can be concluded that a decrease of outer membrane permeability could be possibly due to the reduced expression of porin. The differences in type and content of the OMPs were related to expression levels of the genes encoded by outer membrane efflux protein and porin. Furthermore, it was uncovered that P. aeruginosa resistant to tachyplesin I was associated with reduced outer membrane permeability. It is also noteworthy that the lower expression of proteins should be further identified by MALDI-TOF-TOF. The porin-deficient strains in the identified OMPs might be generated to explore their resistance to tachyplesin I.

The drug resistance mechanism of biofilm formation plays a pivotal role in P. aeruginosa. In a previous study, PA-88 formed an easier biofilm than the original strain PA27853, and the content of extracellular polysaccharide was also higher than that of the original strain.15 The similar results were also found in the PA-99 and PA1.2620 (data were not shown). When the content of exopolysaccharide reaches a certain concentration, it can promote the formation of biofilm. In the present study, we found that a number of genes, formerly known to be associated with biofilm formation, in tachyplesin-resistant PA-99 strain were differentially expressed. The chemosensory pili system protein (chpC) gene was upregulated. The biofilm formation pathway involves the quorum sensing, bacterial secretion system, chemosensory system, etc. The results indicated that P. aeruginosa resistant to tachyplesin I could be involved in the expression of biofilm-related genes. The global interactions of the proteins encoded by the biofilm-related genes were further analyzed.

Some bacteria evade host defense peptides by employing protease-mediated degradation. The PA-99 and PA-88 exhibited the increased levels of extracellular proteolytic activity and reduced the antimicrobial activity of tachyplesin I in previous studies.14,15 Some DEGs were found to encode protease and peptidase in the present study, such as upregulated genes encoding protease (PA3913) and encoding peptidase (PA2658). The upregulated gene PA2315 and downregulated gene PA5514 both encoded beta-lactamase. The downregulated gene PA2764 (log2FC=−16.39) encoded alpha/beta hydrolase superfamily. In the DEGs, some genes that encoded the peptidase S9, S49, and M23 were downregulated. The results revealed that P. aeruginosa resistant to tachyplesin I was associated with the expression of genes, which encoded protease, hydrolase, and peptidase.

To improve the stability of peptides against proteolytic digestion, various chemical modifications, including the use of D-amino acids, cyclization, acetylation, and utilization of delivery systems for AMPs, have been presented.33 Tachyplesin peptides and their cyclized analogs are differentiated by only one or two amino acid residues. A study showed that tachyplesin peptides and cyclized analogs have similarly potent antimicrobial and anticancer properties, while the backbone cyclization improves their stability and therapeutic potential.34

In summary, PA1.2620 strain resistant to tachyplesin I may cause that several genes to be differentially expressed, especially DEGs resistant to drugs. We further verified that P. aeruginosa resistant to tachyplesin I was associated with a decrease of outer membrane permeability, loss of outer membrane porin, and differences in OMPs. The resistance mechanism of P. aeruginosa to tachyplesin I was found to be mainly related to the reduced entrance of tachyplesin I into the bacterial cell through efflux pumps and decrease of outer membrane permeability, which was also similar to antibiotic resistance in P. aeruginosa. The result also indicates that PA-99 strain has many diverse mechanisms including biofilm formation, protease degradation to resist killing by tachyplesin I. However, most changes in the transcriptomes of resistant mutants are the consequence, not the cause of resistance. Therefore, the acquisition of resistance should further be determined the mutations (including gene amplification) involved by whole-genome sequencing to the PA-99 mutant and PA1.2620 strain. Once this is known, transcriptomics may help to track second-order effects as well as to analyze the effect of acquiring resistance on the bacterial physiology.