Genomic characterization of the C. tuberculostearicum species complex, a prominent member of the human skin microbiome

Abstract

Amplicon sequencing data combined with isolate whole genome sequencing have expanded our understanding of Corynebacterium on the skin. Healthy human skin is colonized by a diverse collection of Corynebacterium species, but Corynebacterium tuberculostearicum predominates on many skin sites. Our work supports the emerging idea that C. tuberculostearicum is a species complex encompassing several distinct species. We produced a collection of genomes that help define this complex, including a potentially new species we term Corynebacterium hallux based on a preference for sites on the feet, whole-genome average nucleotide identity, pangenomic analysis, and growth in skin-like media. This isolate collection and high-quality genome resource set the stage for developing engineered strains for both basic and translational clinical studies.

Keywords: Corynebacterium; genomics; pangenome; tuberculostearicum.

Fig 1

Corynebacterium species relative abundance in normal human skin microbiome. (a) Relative abundance of the 15 major Corynebacterium species across 14 skin sites: sebaceous (back, Ba; occiput, Oc; external auditory canal, Ea; retroauricular crease, Ra; manubrium, Mb; glabella, Gb), moist (inguinal crease, Ic; antecubital crease, Ac), dry (hypothenar palm, Hp; volar forearm, Vf), foot (toe nail, Tn; toe web, Tw; plantar heel, Ph), and nares (N). Relative abundances were determined by sequencing of the V1–V3 region of the 16S rRNA gene and subsetting to Corynebacterium reads. (b) Percent of total bacterial reads attributed to Corynebacterium and C. tuberculostearicum in each skin habitat. Of the six ASVs assigned to C. tuberculostearicum, mean relative abundance across skin habitats.

Fig 2

A maximum-likelihood phylogenetic tree of C. tuberculostearicum species complex genomes from this study and publicly available, calculated from 1,315 core gene cluster alignments. Bootstrap values (located along internal nodes) were calculated from 1,000 replicates. Clustering was generated using GET_HOMOLOGUES OrthoMCL v1.4 option with minimum coverage 90% in BLAST pairwise alignments. The tree was rooted on outgroup C. accolens ATCC 49725. On the right of tree, boxes depict site (body site locations defined in Fig. 1) and individual (HV) from which each isolate was cultured. Sites are colored by niche type, with moist in shades of green; feet in shades of orange; dry in pink; sebaceous in lavender; and nares in blue. Individuals are randomly but consistently colored.

Fig 3

The C. tuberculostearicum pangenome. (a) Anvi’o pangenomic map for 28 C. tuberculostearicum genomes (including five NCBI reference genomes). Genomic rings are annotated by skin site and HV metadata and ordered by pyANI ANIb. Genome margins are manually adjusted for clarity. (b) Heap’s Law estimates of pangenome openness for C. tuberculostearicum complex (N = 28) and C. tuberculostearicum species (N = 20) genomes. Rarefaction curves show the total number of genes accumulated with the addition of new genome sequences in random order with 1,000 permutations. Shaded regions represent the 95% CI. A Heap’s law model was fit to the resultant complex and species curves to calculate γ values (0.30 ± 0.01 and 0.25 ± 0.01, respectively). (c) Number of core (belonging to all genomes), accessory (belonging to two or more genomes), and singleton (belonging to only one genome) genes. The expanded pangenome contains 5,451 genes using 90% sequence identity as a cutoff parameter.

Fig 4

C. tuberculostearicum complex pangenome clustering and improved metagenomic read mapping. (a) PCA of orthologous gene clustering. The gene presence/absence data for 25 genomes (including two NCBI references, shown in gray) was analyzed using PCA. Ribotype B genomes are shown as circles; other genomes are triangles. (b) Improvement in shotgun metagenomic read mapping with a 28 member C. tuberculostearicum database as compared to the five member NCBI database. Percent increase in mapped C. tuberculostearicum reads by body site. Each point is an HV. Triangles mark HVs that contributed one or more isolates to the expanded mapping database. (c) Relative abundance of C. tuberculostericum species complex members, including the newly proposed species C. hallux, based on bracken-corrected kraken2 analysis.

Fig 5

Growth phenotypes of select C. tuberculostearicum complex strains in synthetic sweat media. (a) Empirical area under curve comparison of C. tuberculostearicum species complex strains from ribotype A and ribotype B, with biological replicates grouped by color. Strains were grown in BHI + 1% Tween; Sweat media + 0.1% Tween80; Sweat media + 0.1% Tween80 + synthetic lipid mixture. Medium composition is described in further detail in methods. (b–d) Selected growth curves from a representative experiment plotted with standard error. Ribotype B isolates are shown in shades of blue. Ribotype A isolates are shown in shades of red.