Comparative Transcriptome Analysis Reveals Differentially Expressed Genes Related to Antimicrobial Properties of Lysostaphin in Staphylococcus aureus

Comparative transcriptome analysis and de novo short-read assembly of S. aureus Newman strains revealed significant transcriptional changes in response to the exposure to triple-acting staphylolytic peptidoglycan hydrolase (PGH) 1801. Most altered transcriptions were associated with the membrane, cell wall, and related genes, including amidase, peptidase, holin, and phospholipase D/transphosphatidylase. The differential expression of genes obtained from RNA-seq was confirmed by reverse transcription quantitative PCR. Moreover, some of these gene expression changes were consistent with the observed structural perturbations at the DNA and RNA levels. These structural changes in the genes encoding membrane/cell surface proteins and altered gene expressions are the candidates for resistance to these novel antimicrobials. The findings in this study could provide insight into the design of new antimicrobial agents.

In this study, we have identified the genetic mechanisms of resistance to peptidoglycan hydrolases via repeated exposure of S. aureus Newman_2010 strain (cultured wild-type (WT) S. aureus Newman in the year of 2010) to sublethal concentrations of a genetically engineered, triple-acting, staphylolytic, peptidoglycan hydrolase (PGH1801) [2,12]. The resultant mutant strain S. aureus Newman 1801_2010 has resistance to lysostaphin with a >2-fold increase in minimum inhibitory concentration (MIC) relative to the wild-type strain Newman_2010 without exposure to any PGH [2,[13][14][15]. This phenotype does not appear to Antibiotics 2022, 11, 125 2 of 13 be caused by the resistance mechanisms described previously because no changes were found in the known resistance genes (data not shown).
Environmental conditions experienced by multiple generations of S. aureus New-man_2010 cultured with PGH1801 could alter this organism's differential survival and reproduction characteristics. In response to sublethal concentrations of triple-acting PGH1801, S. aureus Newman could produce a unique phenotypic adaptation. This property is known as phenotypic plasticity [16]. Such 'transgenerational' plasticity could provide a competitive advantage in growth resulting in long-term environmental variation [17][18][19][20].
Next-generation sequencing allows the study of any organism at the genomic and transcriptomic levels. Exploring whole-genome expression using next-generation sequencing and RNA-seq provides a more comprehensive understanding than just looking at DNA primary structure, since RNA-seq captures actively transcribed regions and, therefore, can ascertain the molecular basis of a phenotype [21][22][23][24][25]. In this study, we have applied DNAand RNA-seq analyses by comparing genome-wide single nucleotide polymorphism (SNP) and insertion/deletion (InDel) results, as well as by measuring the differential expression levels of thousands of genes simultaneously between the strains of wild-type Newman_2010 and mutant 1801_2010 to identify mutations associated with the above resistant phenotypes. Meanwhile, our comparative transcriptome analysis and de novo short-read assembly revealed the transcriptional changes of S. aureus in response to lysostaphin treatment. All the pooled RNA-seq data were used to compare and establish genome-wide SNP and InDel results between the mutant 1801_2010 and WT Newman_2010. Reverse transcription qPCR was used to confirm the RNA-seq results. Computational analyses, including gene ontology and KEGG pathway enrichment, were also employed. We aim to identify the mutations which provide resistance to the triple-acting PGH. The findings in this work could provide insights into the design of new antimicrobial agents.

Bacterial Strains, Culture Conditions and Mutant Isolation
S. aureus strain Newman_2010 (wild-type) was used in this study [2]. A mutant strain 1801_2010 was developed by repeated exposure of the wild-type strain to sublethal concentrations of the triple-acting PGH1801, a fusion protein containing 3 active domains [2,12,14]. Briefly, S. aureus Newman 2010 bacteria were incubated with two-fold serial dilutions of PGH1801 (10, 5, 2.5, 1.25, and 0.63 µg/well). On day 1 of the MIC assay, bacteria growing in the first turbid well (a sublethal concentration) next to clear wells (lethal concentration) were selected as an inoculate for next MIC assay. Over this time of repeated exposures, bacteria developed more and more resistance to PGH1801. After 10 cycles of exposure of S. aureus Newman_2010 to sublethal concentrations of PGH1801, a mutant strain (1801_2010) with a 2-fold increase in MIC was isolated. Both wild-type and mutant S. aureus strains were grown at 37 • C on tryptic soy agar (TSA) or in tryptic soy broth (TSB). Growth curves of the strains were determined and shown in Supplementary Figure S1. The wild-type and mutant strains displayed a similar growth pattern with the doubling times 47.9 ± 0.431 and 48.1 ± 0.458 min, respectively. using an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) with an RNA 6000 Nano kit. All the samples showed the RNA integrity numbers (RIN) above 8. The efficiency of DNase I treatment of the RNA samples was assessed by PCR amplification of a S. aureus house-keeping gene (gyrA) with a positive (S. aureus genomic DNA) and negative control (nuclease-free water).

RNA-seq Library Construction and Sequencing
Prior to RNA-seq library preparation, ribosomal RNA was depleted from the total RNA using a Ribo-Zero rRNA Removal Kit for Gram-positive Bacteria (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. A Qubit RNA HS Assay and Agilent Bioanalyzer with an RNA 6000 Pico kit were used to assess the quantity and quality of the rRNA-depleted RNA samples. Six libraries (2 strains × 3 biological replicates) were constructed using a TruSeq Stranded mRNA Sample Preparation kit (Illumina, San Diego, CA, USA) following the manufacturer's recommendations. The libraries were sequenced on an Illumina Miseq platform using 300 base-length read chemistry in a paired-end mode.

Reverse Transcription Quantitative PCR (RT-qPCR) Analysis
RT-qPCR was performed on selected genes based on the relative transcription levels (logFC > 2.0) obtained from RNA-seq results. Primers were designed using the NCBI primer designing tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi) (accessed on 10 January 2022). Reverse transcription reactions were carried out using random primers and SuperScript II Reverse Transcriptase (Life Technologies, Carlsbad, CA, USA). Quantification of cDNA was performed on a 7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA) using a SYBR ®® Green PCR Master Mix. The gyrA gene was used as a reference for data normalization. Housekeeping genes gmk and pta were included as controls to ensure data reliability. All the samples were analyzed in three biological and technique triplicates. Relative gene expression levels were computed by the 2 −∆∆CT method, where ∆∆CT = ∆CT (mutant) − ∆CT (WT), ∆CT = CT (target gene) − CT (gyrA), and CT is the threshold cycle value of the amplified gene.

Bioinformatics Analysis
Computational analyses, including gene ontology and KEGG pathway enrichment analysis, were performed [26]. Pooled RNA-seq data were used to compare genome-wide SNP and InDel results between the mutant 1801_2010 and WT Newman_2010 strains. The details are referenced [27].

Data Accession Numbers
Project accessions in NCBI: PRJNA235858 and PRJNA235865 for S. aureus newman_2010 and 1801_2010 strains are available, respectively. DNA sequence accession numbers are SRX478042 (Newman_2010) and SRX478056 (1801_2010). The whole-genome sequence of S. aureus Newman_WT (GenBank accession: NC_009641) was used as a reference for genome mapping.

Integrated Genome and Transcriptome Sequencing for Identification of Genetic Variants
In this study, we applied a Circos diagram ( Figure 1) to illuminate the whole-genome and transcriptome sequencing in order to graphically represent genetic variations and association with their corresponding phenotypes across S. aureus Newman WT and mutant strains. Figure 1 shows several genetic differences between mutant 1801_2010 and WT Newman_2010 compared to the reference genome of Newman_WT (NC_009641). association with their corresponding phenotypes across S. aureus Newman WT and mutant strains. Figure 1 shows several genetic differences between mutant 1801_2010 and WT Newman_2010 compared to the reference genome of Newman_WT (NC_009641).

InDels, Gaps, and SNPs Identified
There was a total of 397 SNPs, 1925 gaps, and 26 InDels detected in mutant 1801_2010 and 360 SNPs, 2185 gaps, and 21 InDels in WT Newman_2010 ( Figure 2) compared to New-man_WT via using the EDGE bioinformatics software [27]. RNA-seq data were used to verify these results. The purpose of this study is to investigate whether these genomic alterations contribute to gene expression. There are six rings in Figure 1: the inner two rings are the genetic difference of InDels, the middle two rings are SNP differences. The second ring from the outside is the list of all genes. All genetic variations of SNPs and InDels from WT Newman_2010 and mutant 1801_2010 were mapped relative to the reference genome sequence of Newman_WT (NC_009641). Several important genes have been identified which could be responsible for the phenotype. Table 1 lists the polymorphic changes in the mutant compared to Newman_WT. Of the 21 significant SNPs, 18 were associated with 4 individual genes (NWMN_0305, NWMN_0306, NWMN_0979, and

InDels, Gaps, and SNPs Identified
There was a total of 397 SNPs, 1925 gaps, and 26 InDels detected in mutant 1801_2010 and 360 SNPs, 2185 gaps, and 21 InDels in WT Newman_2010 ( Figure 2) compared to Newman_WT via using the EDGE bioinformatics software [27]. RNA-seq data were used to verify these results. The purpose of this study is to investigate whether these genomic alterations contribute to gene expression. There are six rings in Figure 1: the inner two rings are the genetic difference of InDels, the middle two rings are SNP differences. The second ring from the outside is the list of all genes. All genetic variations of SNPs and InDels from WT Newman_2010 and mutant 1801_2010 were mapped relative to the reference genome sequence of Newman_WT (NC_009641). Several important genes have been identified which could be responsible for the phenotype. Table 1 lists the polymorphic changes in the mutant compared to Newman_WT. Of the 21 significant SNPs, 18 were associated with 4 individual genes (NWMN_0305, NWMN_0306, NWMN_0979, and NWMN_1288). The remaining three SNPs were distributed in the intergenic regions: two of them were between the lctp and Spa genes and one was located between the recR and tmk genes. Table 2 presents the 15 InDels identified in mutant 1801_2010 compared to Newman_WT. Two genes (NWMN_1308 and NWMN_1410) contained 5 InDels. The remaining 10 InDels were distributed in 5 important intergenic regions: 4 of them were in the intergenic region between the NWMN_0810 and NWMN_0811 genes; 2 were in the intergenic region between the NWMN_2500 and ldh genes; 2 of the InDels were in the intergenic region between the citZ and aapA genes; 2 of them were in the region between dnaA and rpmH. In Tables 1  and 2, two different assembly methods were compared to identify genetic variations: one was based on sequence reads and the other was based on contigs. The results from both methods agreed with each other very well.
the intergenic region between the NWMN_0810 and NWMN_0811 gen intergenic region between the NWMN_2500 and ldh genes; 2 of the In intergenic region between the citZ and aapA genes; 2 of them were in the dnaA and rpmH. In Tables 1 and 2, two different assembly methods w identify genetic variations: one was based on sequence reads and the oth contigs. The results from both methods agreed with each other very well

Identification of Up-Regulated and Down-Regulated Genes
After normalization, the DESeq2 tool [28] was employed for quantitative analysis of RNA-seq data and identification of the genes differentially expressed between the mutant and WT strains. In total, 1091 differentially expressed genes (DEGs) (padj < 0.001) were obtained (Supplementary Table S1): 18 of them were significantly up-regulated (padj < 0.001, log2 FC > 0.7, FDR < 0.05) and 6 genes were significantly down-regulated (padj < 0.001, log2 FC < −0.7, FDR < 0.05). The summary of DEGs between the mutant 1801_2010 and WT Newman_2010 is shown in Table 3. Comparative transcriptome analysis and de novo short-read transcriptome assembly revealed that significant transcriptional changes in response to the triple-acting PGH 1801 were associated with membrane, cell wall, and their related genes (e.g., amidase, peptidase, holin, and phospholipase D/transphosphatidylase).

Function Ontology and KEGG Pathway Enrichment Analyses of DEGs
To annotate the potential functions of the DEGs between WT and mutant strains, DEGs with >2-fold expression change were assigned to different KEGG pathways. All KEGG pathways were analyzed as shown in Figure 3. In these pathways, energy production and conversion, amino acid transport and metabolism, nucleotide transport and metabolism, intracellular trafficking, secretion, and vesicular transport, signal transduction and mechanisms were the most enriched pathways either up-or down-expressed between WT Newman_2010 and mutant 1801_2010 strains.

RT-qPCR Confirmation
To verify the most differentially expressed genes obtained from RNA-seq datasets of WT Newman_2010 and mutant 1801_2010, we performed reverse transcription qPCR experiments and compared these results with RNA-seq. Figure 4 demonstrated that the changes of gene expression from RT-qPCR correlated well with the transcriptome profiling from RNA-seq. Moreover, these results were consistent with the observed structural changes at the DNA and RNA levels (Tables 1 and 2). These structural changes in the genes encoding membrane/cell surface proteins and the perturbation in gene expression are potential candidates responsible for resistance to these novel antimicrobials.

Function Ontology and KEGG Pathway Enrichment Analyses of DEGs
To annotate the potential functions of the DEGs between WT and mutant strains, DEGs with >2-fold expression change were assigned to different KEGG pathways. All KEGG pathways were analyzed as shown in Figure 3. In these pathways, energy production and conversion, amino acid transport and metabolism, nucleotide transport and metabolism, intracellular trafficking, secretion, and vesicular transport, signal transduction and mechanisms were the most enriched pathways either up-or down-expressed between WT Newman_2010 and mutant 1801_2010 strains.

RT-qPCR Confirmation
To verify the most differentially expressed genes obtained from RNA-seq datasets of WT Newman_2010 and mutant 1801_2010, we performed reverse transcription qPCR experiments and compared these results with RNA-seq. Figure 4 demonstrated that the changes of gene expression from RT-qPCR correlated well with the transcriptome profiling from RNA-seq. Moreover, these results were consistent with the observed structural changes at the DNA and RNA levels (Tables 1 and 2). These structural changes in the genes encoding membrane/cell surface proteins and the perturbation in gene expression are potential candidates responsible for resistance to these novel antimicrobials.

Discussion
Our previous genome sequencing work (https://www.ncbi.nlm.nih.gov/sra/?term=SRX478056) (accessed on 12 January 2022) has been trying to identify the genomic changes in Staphylococcus aureus that confer resistance to peptidoglycan hydrolase antimicrobial enzymes [15,29]. In this study, we combined RNA-seq and genome sequence data of WT Newman_2010 and mutant 1801_2010 strains by comparing genome-wide SNP and InDel results. Additionally, we applied RT-qPCR to confirm these results. To reveal the differences between the positive genes and other genes, computational analyses, including gene ontology and KEGG pathway enrichment, were also employed. These differences between mutant 1801_2010, WT Newman_2010, and Newman_WT strains in Figure 1 indicate these genetic changes are potentially responsible for phenotypic variation. However, a phenotypic trait caused by genetic sources of variation could include additive variance, dominant variance, environmental variance (e.g., organismal adaptation), and epistatic variance [30,31].
In Tables 1 and 2, and Figure 2, we observed the genetic variations which could be responsible for the mechanism of resistance to peptidoglycan hydrolase (PGH 1801) in S. aureus. However, these variations are only statistical indicators of a functional effect associated with their genotypic variants because it is uncommon to have a concrete variant

Discussion
Our previous genome sequencing work (https://www.ncbi.nlm.nih.gov/sra/?term= SRX478056) (accessed on 10 January 2022) has been trying to identify the genomic changes in Staphylococcus aureus that confer resistance to peptidoglycan hydrolase antimicrobial enzymes [15,29]. In this study, we combined RNA-seq and genome sequence data of WT Newman_2010 and mutant 1801_2010 strains by comparing genome-wide SNP and InDel results. Additionally, we applied RT-qPCR to confirm these results. To reveal the differences between the positive genes and other genes, computational analyses, including gene ontology and KEGG pathway enrichment, were also employed. These differences between mutant 1801_2010, WT Newman_2010, and Newman_WT strains in Figure 1 indicate these genetic changes are potentially responsible for phenotypic variation. However, a phenotypic trait caused by genetic sources of variation could include additive variance, dominant variance, environmental variance (e.g., organismal adaptation), and epistatic variance [30,31].
In Tables 1 and 2, and Figure 2, we observed the genetic variations which could be responsible for the mechanism of resistance to peptidoglycan hydrolase (PGH 1801) in S. aureus. However, these variations are only statistical indicators of a functional effect associated with their genotypic variants because it is uncommon to have a concrete variant with a "precise" genetic location with measurable statistical "effects". Although we can say that there is a functional effect, in order to determine functionality we must go beyond the identification of the variant phenotype due to a specific SNP-or InDel-based locus. It will be instructive to compare the results of current study with those in the literatures. Mostly, there are genomic loci that could influence the expression level of mRNA and these loci can be physically located close to or far away from the gene that gets regulated. It is not necessary that genetic loci are associated with a SNP or InDel.
Through the Gene Ontology function and KEGG pathway enrichment analyses, we found energy production and conversion, amino acid transport and metabolism, nucleotide transport and metabolism, intracellular trafficking, secretion, and vesicular transport, signal transduction and mechanisms were the most enriched pathways either up-or downexpressed between the wild-type and mutant.
Comparative transcriptome analysis and de novo short-read assembly in this study revealed that the genes with significant transcriptional changes in response to exposure to the triple-acting fusion protein are associated with membrane and cell wall (e.g., amidase, peptidase, holin, and phospholipase D/transphosphatidylase). These results are consistent with the observed nucleotide changes at the DNA level. The nucleotide changes in the genes encoding membrane/cell surface proteins and the alteration of gene expression may contribute to the increased resistance of S. aureus to PGHs. The findings of this study could provide insights into the target genes responsible for PGH resistance and lead to the design of new antimicrobial agents.