Whole-Genome Sequencing of Staphylococcus cohnii Isolated from Healthy Human Skin: Insights into Genomic Features and Antibacterial Potential

1. Introduction

Finding effective antibiotics to treat bacterial infections is a major concern, as antibiotic resistance (AMR) is becoming increasingly common, threatening treatment capabilities. The misuse of antibiotics has led to the rise of multidrug-resistant (MDR) bacteria, reducing the effectiveness of conventional antimicrobial drugs. The World Health Organization (WHO) has warned that the situation of antibiotic resistance is worsening, while the Centers for Disease Control and Prevention (CDC) has declared that humanity is now living in a "post-antibiotic era" [1].

The human skin microbiome contains hundreds of different types of bacteria that serve as the body's first defense against the outside environment [2]. Coagulase-Negative Staphylococci (CoNS) are bacteria generally considered non-pathogenic and are commonly found as part of the human skin microbiome [3]. These bacteria create heat-stable, ribosomal synthesized antimicrobial peptides called bacteriocins, which can have both broad and narrow range inhibitory spectra. These bacteriocin have no antagonistic effects on the cells that make them [4,5].

Staphylococcins produced by Staphylococci display remarkable functional diversity, being effective not only against other staphylococci but also against a wide range of microorganisms. Likewise, "bacteriocin-like inhibitory substances" (BLIS) refer to biomolecules which inhibit widely range of bacteria including both Gram-positive and Gram-negative bacteria [6]. Gram-positive bacteriocins are categorized into three main classes based on their molecular structure and post-translational modifications. Class I, also called lantibiotics, includes molecules under 5 kDa that are thermostable and feature significant post-translational modifications, with distinct amino acid compositions like lanthionine and methyl-lanthionine forming intramolecular rings; this class is further divided into Class Ia, which are polar with a positive net charge, and Class Ib, which are neutral or negatively charged. Class II comprises small molecules (<10 kDa) with minimal or no post-translational modifications and is subdivided into four subclasses based on their structure and mode of action: IIa, pediocin-like peptides with conserved N-terminal YGNG(V/L)XC sequence; IIb, unmodified two-component bacteriocins; IIc, circular peptides; and IId, other single-peptide unmodified bacteriocins [7]. Class III encompasses large proteinaceous bacteriocins (>30 kDa) with complex, thermolabile structures, some of which can also be produced by Gram-negative bacteria [8]. In our previous study, phenotypic assays revealed significant antibacterial activity of the Staphylococci strain habitated from human epidermis against several bacterial species, including Staphylococcus aureus ATCC 29213, Staphylococcus aureus ATCC 43300, Enterococcus faecalis ATCC 29212, Enterococcus faecalis ATCC 51299, Enterococcus casseliflavus ATCC 700327, and Staphylococcus saprophyticus ATCC BAA-750 (data not shown). Thus, the present study aims to: (1) characterize the complete genome of a Staphylococci strain isolated from healthy human skin; (2) identify and analyze the propebacteriocin genes within its genome, with an emphasis on their potential antimicrobial properties; (3) investigate the phylogenetic relationships between the Staphylococci strain and other related Staphylococcus species to elucidate its evolutionary dynamics.

2. Materials and Methods

2.1. Bacterial Strain and Culture Condition

Tested strains were isolated from healthy human skin and cultured under aerobic conditions at 37^oC in BHA medium (Brain heart agar) for 24h in an incubator. Isolated colonies were subjected to Gram staining to differentiate between cocci, bacilli, and fungi, and underwent catalase, tube coagulase tests, and oxidase testing to distinguish between coagulase-positive Staphylococcus (CoPS), coagulase-negative Staphylococcus (CoNS), and Micrococcus spp [9]. Participation was voluntary, and every participant signed a written informed consent form. The study was approved by the Ethics Committee in Biomedical Research, Hue University of Medicine and Pharmacy, Hue University (Approval Code H2023/062; May 20, 2023).

2.1. Genome Sequencing and Assembly

DNA extraction:

Genomic DNA extraction was carried out using the Monarch HMW DNA extraction kit for tissue (NEB, USA). DNA concentration was determined by fluorescence method by Qubit DNA assay (Thermo Fisher Scientific) while the OD260 and OD280 values were determined by Nanodrop spectrophotometer. DNA size was assessed with 1% agarose gel electrophoresis. Extracted genomes must had a concentration of at least 50.00 ng/μL, OD260/OD280 ≥ 1.50. DNA size < 1000 bp will be considered carefully.

Library preparation and sequencing

To overcome the drawbacks of both short- and long-read sequencing, we employed a hybrid assembly strategy that integrates the high accuracy of Illumina reads with the long-range continuity provided by Oxford Nanopore. Paired-end short-read libraries with 150 bp reads were prepared using the NEBNext dsDNA Fragmentase, NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB, USA) while 8000 bp long-read libraries were carried out by Ligation Sequencing Kit (Oxford Nanopore Technologies, USA) according to the manufacturer's instructions. Short-read libraries were sequenced using next-generation sequencing on Novaseq (Illumina) and long-read libraries were sequenced using next-generation sequencing on Flongle (Oxford Nanopore Technologies). The size of short-read and long-read libraries were obtained within 0.4 Gb and 0.15 Gb, respectively.

Genome assembly

To ensure high-quality data, raw sequencing reads were first cleaned and filtered. Short reads were processed with fastp v0.23.1 [10], while long reads were handled using filtlong v0.2.1. Low quality nucleotides and ambiguous bases ("N"s) were removed based on Phred quality scores. Long reads were further filtered by length. The resulting high-quality short and long reads were then assembled into a genome using two independent de novo approaches: Unicycler v0.4.8 [11] and a pipeline consisting of Flye v2.9 – medaka v1.4.3 -polca v4.0.9 [12,13]. The best assembly, chosen by comparing the outputs of these methods, was then selected for further analysis. Assembly quality was assessed using several methods, including Quast v5.2.0 [14] for quantitative metrics and checkM v1.2.1 for genome completeness [15]. To refine the assembly, Dnaapler was used to correct the chromosome start position.

2.4. Gene Content Analysis

Genome annotation

The assembled genome was annotated using Prokka v1.14.6 [16] in genus-specific mode to identify genes and other genomic features. To access the completeness of the genome assembly and gene annotation, BUSCO v5.8.2 [17] was used with “genome” and “proteins” modes and the lineage dataset of “bacteria_odb12”.

Amino acid sequences of the annotated genes were used as input for the COGclassifier (https://github.com/moshi4/COGclassifier) and GhostKOALA [18] software packages to classify them into different functional groups. GhostKOALA classification results were visualized using “ggplot2” and “ggpubr” R packages [19,20].

Bacteriocin encoding gene identification and analysis

BLASTp [21] was used to identify the suspected bacteriocin with custom database of reference bacteriocin sequences from previous studies [22]. The E-value cutoff for BLASTp was set as 0.01. The amino acid sequences having the best BLAST hits based on bit scores, E-value, and nucleotide were considered as hypothetical bacteriocins and used for further analysis.

2.5. Phylogenetic Analysis

16S rDNA phylogenetic tree construction

The 16S rDNA region of the collected isolate was amplified by using universal primer pairs of 27F (5´-AGAGTTTGATCCTGGCTCAG-3´) and 1492R (5´-GGTTACCTTGTTACGACTT-3´) on genomic DNA template as described by Ha et al [23]. Then, the 16S rDNA sequences of the isolate and other Staphylococcus species, including S. haemolyticus, S. epidermidis, S. aureus, S. saprophyticus, S. capitis and S. hominis were aligned using MAFFT v7.526 [24] with the alignment strategy set as “auto”. The alignment file was then used as input for IQ-TREE v2.3.6 (-m MFP -B 5000 -bnni -alrt 1000 -seed 31081997 -st DNA) [25].

Bacteriocin gene phylogenetic construction

To investigate the origins of hypothetical bacteriocin genes, Orthofinder v.3.0.1b1 [26] was used to find orthologs of these genes in 14 other Staphylococcus species, with default parameters. The Staphylococcus species from which 14 genomes were downloaded from EsemblBacteria (https://bacteria.ensembl.org/index.html) include S. agnetis (GCA_003040835), S. argenteus (GCA_003967115), S. arlettae (GCA_900457375), S. aureus (GCA_000638495), S. auricularis (GCA_002902455), S. capitis (GCA_001179785), S. caprae (GCA_000160215), S. carnosus (GCA_003970565), S. chromogenes (GCA_003040195), S. condimenti (GCA_001618885), S. croceilyticus (GCA_002902575), S. debuckii (GCA_003718735), S. delphini (GCA_002369735), and S. devriesei (GCA_002902625). The ortho-groups containing a single gene for each species were concatenated and aligned using MAFFT v7.526 [27] with the alignment strategy set as “auto”. The alignment file was then used as input for IQ-TREE v2.3.6 (-m MFP -B 5000 -bnni -alrt 1000 -seed 31081997 -st AA) [25].

The GFF3 and genome FASTA files of the Staphylococcus species having bacteriocin orthologs that are the most identical to the ones in our isolate were downloaded from EnsemblBacteria. Annotation features of the bacteriocin orthologs were then extracted from the GFF3 file. The extracted annotation features for bacteriocin orthologs and the genome sequences were used to generate gbk files using the “seqret” tools of EMBOSS v.6.6 [28]. The identity of bacteriocin orthologs were visualized using Clinker v0.0.31 [29] and the script “Genbank_slicer.py” (https://github.com/jrjhealey/bioinfo-tools/blob/master/Genbank_slicer.py).

Phylogenetic tree visualization

All generated Maximum-likelihood trees were visualized using the “ggtree” package from the R programming language [30].

3. Results

3.1. 16S Phylogenetic Analysis

Presumptive coagulase-negative Staphylococci isolates were identified from the skin swabs based on exhibiting coagulase-negative and oxidase-negative phenotypes. To identify the Staphylococci isolate (denoted by 148-XN2B 18.2), the 16S rDNA of the strain were amplified and compared with the 16S rDNA sequence of S. haemolyticus, S. epidermidis, S. aureus, S. saprophyticus, S. capitis and S. hominis. The resulting phylogeny tree (Figure 1) showed that strain 148-XN2B 18.2 has a high similarity to Staphylococcus cohnii ATCC 29974 (OR801636) and S. cohnii subsp. cohnii (D83361.1), with a sequence identity of 100%. Additionally, strain 148-XN2B 18.2 also exhibited 100% identity with strains isolated in other studies, including S. cohnii GH137 and S. cohnii GS29, and 99.93% similarity with strain S. cohnii CTB46.

3.2. Genome Analysis of S. Cohnii 148-XN2B 18.2

Regarding short-read sequencing metrics, the total raw reads reached 3,376,387, with an average GC content of 32.7% and the rate of reads with Q30 was 96.5%. The GC content of the genome was 32.7%. For long-read, a total of 47,648 raw reads was achieved, with the rate of reads reaching Q7 was 100%. Genome assembly was performed using both short and long reads, and revealed a total of 4 contigs. The complete genome of S. cohnii 148-XN2B 18.2 was identified to have a single circular chromosome with a total length of 2,768,657 bp and no plasmids (Figure 2). Within the theoretical coding regions, 2,756 genes were predicted, including 2,676 CDS (protein-coding regions), which accounted for the majority of the genome's total length (Table 1). The complete genome sequence has been deposited in NCBI with the project accession SAMN54131608.

3.3. Gene Annotation and Functional Annotation

The functional classification of annotated genes was performed using GhostKOALA, and the results were visualized in Figure 3. A total of 19 functional groups were identified, in which the largest group was Genetic Information Processing (A), comprising over 400 genes, followed by Signaling and Cellular Processes (B) with approximately 200 genes and Carbohydrate Metabolism (C) with over 150 genes. In addition, genes associated with Metabolism (D), Environmental Information Processing (E) and Amino Acid Metabolism (F) were also significantly represented, each containing around 100 genes. Meanwhile, other categories had fewer genes. Notably, the Unclassified group (I) accounted for a notable number of genes, approximately 100.

3.4. Identification of Bacteriocin Genes in S. cohnii 148-XN2B 18.2

A previous study on genome analysis of bacteriocin-producing Staphylococci showed that the their genomes contain various gene encoding for bacteriocin class I and class II [31]. Thus, the genome of S. cohnii were searched for the potential genes encode for bacteriocin class I and class II. The results present in Table 2 indicated 5 potential bacteriocin-related genes identified in the genome of S. cohnii, along with their reference bacteriocin IDs and associated information. Among these, S. cohnii carries one gene for class I and four genes for class II.

3.5. Phylogenetic Relationship and Comparative Genomic Analysis of S. cohnii and S. arlettae

To find the genetic context of S. cohnii isolate, an ortholog phylogenetic tree was constructed between different Staphylococci. The phylogenetic analysis based on the orthogroups of different Staphylococci demonstrates a close evolutionary relationship between tested S. cohnii and S. arlettae, supported by a bootstrap value of 95.2%. This analysis highlights the strong genetic similarities between these two species (Figure 4), which would be investigated for the diversity of their bacteriocin contents.

The comparative genomic analysis of S. cohnii and S. arlettae’ revealed three bacteriocin-related gene clusters with conserved and homologous regions between the two species. The identified clusters ranged in size from 0.3 kb to 2.5 kb, with sequence identity values represented in grayscale for comparative visualization. Colicin V-like gene cluster was identified in S. cohnii (coordinates: 1,734,743–1,735,624) and exhibited high sequence identity with the corresponding region in S. arlettae (coordinates: 1,821,949–1,822,940). Besides, Enterocin B-like gene cluster was observed in S. cohnii (coordinates: 2,216,019–2,215,730) and its homologous region in S. arlettae (coordinates: 2,216,590–2,190,969), which also encoded a Pyocin AP41-like gene. A Bacteriocin 28b-like gene cluster was also identified in S. cohnii (coordinates: 2,386,002–2,386,520) and S. arlettae (coordinates: 2,448,002–2,448,651), with strong sequence similarity (Figure 5).

To provide a broader context for the genomic data of an isolated Staphylococci from healthy human skin, a 16S rRNA phylogenetic tree was constructed between S. cohnii and related Staphylococcus spp, which showed the 148-XN2B 18.2 isolate to be closely related to commensal coagulase-negative S. cohnii strains, such as S. cohnii ATCC29974 and GS29, while differed from other common coagulase-positive Staphylococcus species like S. aureus and S. epidermidis . The functional annotation of S. cohnii 148-XN2B 18.2 genes reveals a genome structured to prioritize essential processes, particularly genetic information processing pathways, signaling and cellular processes, and different metabolism pathways. This focus aligns with the needs of a commensal organism, where maintaining transcriptional, translational, and replicative fidelity is crucial for survival under diverse environmental conditions [32,33,34,35,36,37]. Class II bacteriocins were also identified, comprising Acidocin 8912, Enterocin B, Bacteriocin 28b and Colicin V, which exert antimicrobial effects through mechanisms such as pore formation and membrane destabilization in target bacteria. Identified class II bacteriocins had demonstrated a wide range of target, from Gram-positive bacteria (Acidocin 8912, Enterocin B) [4,38], Gram-negative bacteria (Colicin V) [39], to both targets (Bacteriocin 28b) [40]. The coexistence of these three potent bacteriocin-encoding genes in S. cohnii underscores its ecological adaptability and biotechnological potential in applications such as food preservation and therapeutic development [41]. In addition, a Class I Nif11-like peptide family was also identified, which are notable for their structural stability and potential applications in synthetic biology due to their post-translational modifications [42]. The discovery of such genes highlights S. cohnii’s potential to produce highly stable and functionally diverse antimicrobial peptides.

Based on the orthogroup phylogenetic tree between S. cohnii 148-XN2B and different Staphylococcus spp., a close phylogenetic relationship between S. cohnii and S. arlettae was unveiled. This clustering aligns with previous studies that have identified S. cohnii and S. arlettae as part of a group of coagulase-negative staphylococci (CoNS) frequently isolated from similar host environments, such as human and animal skin microbiomes [43,44]. The comparative genomic analysis emphasized the conservation of the Colicin V-like gene cluster, Enterocin B-like cluster, Pyocin AP41-like gene and Bacteriocin 28b-like cluster between S. cohnii and S. arlettae, which may reflect their shared ecological niches and evolutionary pressures, suggesting that both S. cohnii and S. arlettae may rely on similar mechanisms to inhibit competing bacteria in their environments [45]. Additionally, the functional conservation of these bacteriocin gene clusters may have implications for the development of novel antimicrobial agents. The shared genetic architecture between S. cohnii and S. arlettae highlights their potential as reservoirs for diverse bacteriocins with therapeutic applications.

However, there were some limitations in this study. While the classification tools COGclassifier and GhostKOALA provided valuable insights, they rely on existing reference datasets, which might limit the detection of novel functional genes. Additionally, targeted analysis of underrepresented categories, such as secondary metabolite biosynthesis, could reveal further insights into the organism's biosynthetic potential. Besides, the functional roles of several bacteriocin genes identified in this study have not been experimentally validated, which makes it challenging to confirm their precise contribution to the observed antimicrobial activity. Moreover, translating these findings into clinical applications may face hurdles such as scaling up bacteriocin production, navigating regulatory requirements, and addressing variability in efficacy due to genetic differences among strains or environmental factors. Despite these challenges, future studies focusing on functional assays, optimization of production processes, and exploration of synergistic applications with other antimicrobial agents could significantly advance the utility of S. cohnii-derived bacteriocins.

4. Conclusions

This study provides a comprehensive genomic and functional analysis of S. cohnii 148-XN2B 18.2, emphasizing its potential as a reservoir of bacteriocins. The evolutionary conservation and demonstrated antimicrobial activity of its bacteriocins make it a promising candidate for future applications in addressing antimicrobial resistance.