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

Abstract

Corynebacterium is a predominant genus in the skin microbiome, yet its genetic diversity on skin is incompletely characterized and lacks a comprehensive set of reference genomes. Our work aims to investigate the distribution of Corynebacterium species on the skin, as well as to expand the existing genome reference catalog to enable more complete characterization of skin metagenomes. We used V1-V3 16S rRNA gene sequencing data from 14 body sites of 23 healthy volunteers to characterize Corynebacterium diversity and distribution across healthy human skin. Corynebacterium tuberculostearicum is the predominant species found on human skin and we identified two distinct C. tuberculostearicum ribotypes (A & B) that can be distinguished by variation in the 16S rRNA V1-V3 sequence. One is distributed across all body sites and the other found primarily on the feet. We performed whole genome sequencing of 40 C. tuberculostearicum isolates cultured from the skin of five healthy individuals across seven skin sites. We generated five closed genomes of diverse C. tuberculostearicum which revealed that C. tuberculostearicum isolates are largely syntenic and carry a diversity of methylation patterns, plasmids and CRISPR/Cas systems. The pangenome of C. tuberculostearicum is open with a core genome size of 1806 genes and a pangenome size of 5451 total genes. This expanded pangenome enabled the mapping of 24% more C. tuberculostearicum reads from shotgun metagenomic datasets derived from skin body sites. Finally, while the genomes from this study all fall within a C. tuberculostearicum species complex, the ribotype B isolates may constitute a new species.

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 (N)ares. Relative abundances 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 1315 core gene cluster alignments. Bootstrap values (located along internal nodes) were calculated from 1000 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 Figure 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 Corynebacterium tuberculostearicum pangenome. (A) Anvi’o pangenomic map for 28 C. tuberculostearicum genomes (including 5 NCBI reference genomes). Genomic rings are annotated by skin site and HV (healthy volunteer) metadata and ordered by pyANI average nucleotide identity (ANIb). Genome margins are manually adjusted for clarity. (B) Heap’s Law estimate of pangenome openness for 28 genomes. A rarefaction curve showing the total number of genes accumulated with the addition of new genome sequences in random order with 1000 permutations. Shaded regions represent the 95% confidence interval. A Heap’s law model was fit to the resultant curve to calculate k and γ values (1977±38.0 and 0.30±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 5451 genes using 90% sequence identity as a cutoff parameter.

FIG 4

C. tuberculostearicum complex pangenome clustering and improved metagenomic read mapping. (A) Principal components analysis of orthologous gene clustering. The gene presence/absence data for 25 genomes (including two NCBI references, shown in gray) was analyzed using principal components analysis. 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 5 member NCBI database. Percent increase in mapped C. tuberculostearicum reads by body site. Each point is a healthy volunteer. Triangles mark healthy volunteers that contributed one or more isolates to the expanded mapping database.

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 Brain Heart Infusion (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.