Intraspecies warfare restricts strain coexistence in human skin microbiomes

Abstract

Determining why only a fraction of encountered or applied strains engraft in a given person's microbiome is crucial for understanding and engineering these communities. Previous work has established that metabolic competition between bacteria can restrict colonization success in vivo, but other mechanisms may also prevent successful engraftment. Here we combine genomic analysis and high-throughput agar competition assays to demonstrate that intraspecies warfare presents a significant barrier to strain coexistence in the human skin microbiome by profiling 14,884 pairwise interactions between Staphylococcus epidermidis isolates cultured from 18 people from 6 families. We find that intraspecies antagonisms are abundant, mechanistically diverse, independent of strain relatedness and consistent with rapid evolution via horizontal gene transfer. Critically, these antagonisms are significantly depleted among strains residing on the same person relative to random assemblages, indicating a significant in vivo role. Wide variation in antimicrobial production and resistance suggests trade-offs between these factors and other fitness determinants. Together, our results emphasize that accounting for intraspecies warfare may be essential to the design of long-lasting probiotic therapeutics.

Extended Data Figure 1:. Antimicrobial BGCs and defense genes are enriched in accessory genomes.

(A) We constructed a pangenome for S. epidermidis isolates in our collection and annotated genes and biosynthetic gene clusters (BGCs). Antagonism-related BGCs were more frequent in accessory genomes, whereas metabolism-related BGCs were found in the core genome. (B) We identified regions of the genome with copy number variation within a lineage, indicating recently gained or lost regions. We compared the relative enrichment of COG terms and hits from functional annotation databases among genes in the core genome and recent gain/loss regions. (C) As in B, but comparing genes in the core and accessory genomes. COG categories denoted in blue are metabolism associated, whereas antimicrobial production and resistance genes may be found in red categories. Significant enrichment in a two-sided Fisher’s Exact test is indicated by an asterisk (one for naive P-value < 0.05, two for Bonferroni corrected P-value < 0.05), and the number of genes annotated in each category is indicated to the right of the figure.

Extended Data Figure 2:. Representative isolates capture the S. epidermidis lineage diversity of individuals and families.

In a previous study, up to 96 isolates of S. epidermidis were cultured from each of 23 subjects at each timepoint sampled during a 1.5-year time period. Six families were selected for inclusion in this study. Relative lineage abundance, as determined from isolates, is plotted for each sample, organized by subject and family. PHLAME was applied to metagenomic samples from the same subject (denoted as MG), which confirmed that lineage abundances estimated by culturing did not have significant culture bias; note that uncalled portions of metagenomic samples (space above bars) can reflect uncertainty in PHLAME due low coverage rather than the true absence of lineages. We used isolate-based lineage abundances because deep metagenomic coverage of S. epidermidis was not available for every sample. Sample labels indicate the family, subject, and sampling timepoint (e.g. 1AA3 indicates Family 1, subject 1AA, timepoint 3). The number of S. epidermidis colonies isolated from each person is indicated above each bar. Lineages represented in the TSA antagonism screen in this study are indicated in color, while unrepresented lineages are indicated in greyscale. Some additional lineages were included in a follow-up screen M9+S screen (see Methods); these are indicated by darker grey colors (denoted by *). Colors do not carry over across families, lineages are generally restricted to a family. Across all samples, the average proportion of lineage composition represented by isolates included in the TSA antagonism screen was 96%. Subject 8PB reported antibiotic use between timepoints, which may explain the difference between timepoints.

Extended Data Figure 3:. Antagonism and sensitivity generally cluster at the lineage level.

We measured 21,025 pairwise interactions between 145 isolates that passed inclusion criteria (15 replicate cultures were also included, totaling 25,600 interactions). The heatmap depicts the AUC of the ZOI calculated for pairwise interactions between S. epidermidis and other skin microbiome isolates in TSA. Each row and column indicate an isolate (in contrast to Fig. 2a where rows indicate lineages). Multiple isolates from a lineage are denoted by decimals (e.g. 20.1, 20.2, 20.3). Isolates are sorted by phylogenetic similarity, with dashed lines to indicate species boundaries. Each isolate of species other than S. epidermidis (labels > 200, see x-axis for details) was treated as a separate lineage. Interactions were shared among all isolates of the same lineage in ~97% of cases, so S. epidermidis isolates were dereplicated into lineages for other analyses unless otherwise noted. Nine isolates were added after the preliminary TSA screen, and thus are lacking some data points, indicated by brown squares. Unless otherwise noted, these isolates with partial data were excluded from statistical analysis.

Extended Data Figure 4:. Antagonism frequency and strength show limited variation at different phylogenetic distances.

Rather than depicting genomic relatedness using a tree (Fig. 2a), we calculated pairwise distances by different sequence similarities to understand whether antagonism and relatedness are correlated. (A) We binned interactions between pairs of lineages with similar average amino acid identities, then plotted the frequency of antagonism among pairs in each bin. (B) As in A, but binned according to average nucleotide identity. (C) As in B, but only for lineages of S. epidermidis, not other species. Antagonism frequency does not significantly vary with genomic distance in our collection, as both R2 and P are weak. (D) We plotted the strength of antagonism (AUC of ZOI) for each antagonistic interaction (excluding non-antagonistic interactions) against the average amino acid identity of the pair. (E) As in D, but against the average nucleotide identity. (F) As in E, but only for lineages of S. epidermidis, not other species. There is a significant but small increase in antagonism strength at high nucleotide identities among S. epidermidis. All P-values and R2 values report the fit of a linear model evaluated by a two-sided ANOVA. (G) The frequency of antagonism does not significantly differ between members of the same or different S. epidermidis phylogroup, according to a two-sided permutation test with shuffled group labels. The phylogeny at right depicts S. epidermidis phylogroup boundaries. Taken together, these results indicate that antagonism occurs at all levels of genomic relatedness, even between closely related lineages.

Extended Data Figure 5:. Warfare by S. epidermidis strains is mediated by diverse and unannotated effectors.

In order to check for shared mechanisms, the binary antagonism matrix was clustered by constructing dendrograms for antagonism and sensitivity (Methods). White boxes indicate antagonism by the producer spot against the receiver lawn. Red arrows indicate the lineages in our collection chosen for our mechanism screen, which include 18 lineages with high antagonism frequencies plus 2 negative controls (which lacked antimicrobial production genes but contained suspected prophage, denoted by yellow arrows). We chose representative receiver bait lineages using the dendrogram of sensitivity to ensure that all interaction patterns were captured. Note that for this analysis, unlike other analyses in this study (Methods), antagonisms from idiosyncratic isolates that exhibit intra-lineage variation were not excluded. Lineages denoted by cyan arrows did not inhibit many S. epidermidis, with the exception of idiosyncratic lineages 20, 37, and 58, which cluster with some S. hominis and S. capitis isolates, denoted by a cyan line. It is possible that the antimicrobials behind these interactions are intended to target non-epidermidis species. (B) We screened cells and filtered supernatants against a representative array of bait lawns chosen based on interaction profiles from above. Filtrates treated with heat or proteinase to degrade proteins and peptides, desalting to remove small molecules, or 1:10 dilution were also tested. Filtrates exhibited varied patterns, including antagonisms likely mediated by small peptide effectors (isolate 30.1) and some cell-associated antagonisms induced by neighbors (isolate 34.2) or both. (C) Some cell-associated antagonisms are iron-mediated, as indicated by the elimination of the interaction in media supplemented with excess FeCl3 (see isolate 51.1). All tested isolates displaying iron-suppressible antagonisms cells contained a biosynthetic gene cluster (BGC) for a product similar to pulcherriminic acid, an iron-sequestering pigment molecule which sabotages the iron supply. (D) We annotated putative bacteriocin BGCs and found BGCs for known compounds (including gallidermin-like, epicidin-like, and lactococcin-like BGCs) as well as unknown compounds. However, we did not find a clear correspondence between BGC presence and the phenotype classified in our mechanism screen, potentially due to multiple inducible effectors being encoded in the genome, incomplete expression of BGCs in different media or environments, or targets not represented in our screen (e.g. Cutibacterium acnes). One isolate, 32.1, had an active prophage that occasionally produced distinct plaques in this assay, and was later isolated using mitomycin-C induction.

Extended Data Figure 6:. Antagonism differences between TSA and M9+S do not affect depletion of on-person antagonism.

(A) We overlaid interactions from TSA onto the S. epidermidis interaction matrix generated in M9+S media, using color to indicate which media the antagonism appeared in. In total, 91% of interactions are the same in both media conditions, though antagonisms did vary between the two media, such that 23% of antagonisms observed in TSA were also observed in M9+S, and 61% of antagonisms observed in M9+S were observed in TSA. (B-D) Despite these differences, we observed a significant depletion of on-person antagonism (ΔAF, see Fig. 3c, Methods) in permutation tests that compare our observed data compared to simulated data with shuffled lineage compositions on each subject. Significance P-values are derived from two-sided permutation tests with shuffled lineage labels (Methods), without multiple hypothesis correction.

Extended Data Figure 7:. Antagonistic interactions are depleted in samples, individuals, and families.

(A) Antagonism frequencies (AF) between co-resident lineages in samples appears depleted relative to the frequency across all S. epidermidis lineages (red). When lineages were present in a sample but were absent from our screen, we calculated the range of AF that would be expected based on the average rate of antagonism between any two isolates. Error bars depict the upper bound of the AF range in a given sample. This is likely an overestimate of the range, since on-person antagonism frequencies are lower than population-wide antagonism frequency. (B) We took the mean AF across samples of the same subject, and still found that mean within sample AF appeared to be depleted relative to the expected AF. Reproduced from Fig. 3a. (C-F) Since t-tests are invalid for non-independent samples, we performed permutation testing by simulating scenarios that retained the structure in the interaction matrix and the composition of samples (e.g. uneven lineage distribution). For each simulation, we shuffled lineage labels in the relative abundance matrix and the p-value represents the fraction of simulations with less on-person antagonism than observed (black arrow indicates observed value of ΔAF). (C) The rows of the composition matrix were shuffled once across all lineages in our collection. This simulation tests for depletion of antagonism while maintaining the same degree of lineage sharing seen between individuals in the same family. Reproduced from Fig. 3c. (D) The rows of the composition matrix were shuffled, AF values were calculated for samples of the first subject, then rows were shuffled again before each new subject. These simulations broke the relatedness between individuals in the same family, while retaining similarity between samples from the same subject, thereby testing depletion of antagonism under the assumption that lineages are acquired from the environment independent of family. (E) To test for depletion of antagonism under the assumption that one can only acquire strains from a family member, shuffling was only performed on submatrices composed only of lineages on a family were shuffled, by row, as in C. (F) As in F, but without maintenance of lineage sharing between family members; submatrices shuffled but by once per subject. Note that less shuffled conditions (E,F) will definitionally have higher p-values because of fewer possible permutations. Under all tested assumptions, antagonism is depleted among strains present on the same person at the same time. Significance P-values are derived from two-sided permutation tests with shuffled lineage labels (Methods), without multiple hypothesis correction.

Extended Data Figure 8:. Subjects colonized by antagonists have different lineage compositions than family members.

Relative abundance plots where each lineage is colored by the proportion of other lineages it antagonizes on TSA, with blues representing non-antagonistic lineages and reds representing highly antagonistic lineages. Lineages are labelled by number in white text and are plotted in the same order as in Extended Data Fig. 2. Grey lineages have missing data, and therefore an antagonism proportion cannot be calculated. Note that antagonism proportion is relative to other lineages across the study, not just within a family. Some families contain few antagonistic lineages (e.g. family 1), whereas most families contain at least one highly antagonistic lineage.

Extended Data Figure 9:. Phenotypic variation among S. epidermidis.

(A) Lineages were grouped into phylogroups according to phylogenetic structure as described previously, see Extended Data Fig. 4g. On many subjects, multiple phylogroups coexist on the face at the same time. (B) The agrD sequence was annotated for each lineage on each subject and classified into quorum sensing types. Again, we found that on many subjects, multiple agr types coexist on the face at the same time. Co-residence between different phylogroups and between incompatible quorum sensing types suggest that there is not competitive exclusion on the basis of phylogenic dissimilarity or quorum sensing. (C) We measured growth rate in liquid culture in TSB for all isolates in our collection (Methods). We calculated the median growth rate for each isolate in a lineage. Bars denote standard deviation (from technical duplicates of biological triplicates). (D) Lineages that grow faster than other lineages in vitro (green) reach higher abundance on individuals (blue), Spearman’s rank correlation. Note that growth in vivo may vary significantly due to differences between nutrient composition of skin and TSB. (E) Lineages that antagonize a higher fraction of other lineages (red) do not grow faster or slower than other lineages (green), Spearman’s rank correlation. (F) The growth rate for each isolate of idiosyncratic lineages 20, 37, and 58 (technical duplicates of biological quadruplicates) is shown. Growth rate varies between vraFG mutants (indicated with *) and wild-type isolates of each lineage, but without a directional trend. Reversion mutants of 37 have faster growth rates than 37.3, leaving the basis of the vraFG mutations unresolved.

Extended Data Figure 10:. Rare cases of derived sensitivity to antimicrobials and phage resistance.

Four “Idiosyncratic” isolates (20.3, 37.3, 49.5, and 70.3) exhibited large differences in antimicrobial sensitivity compared to other members of their respective lineages. Two idiosyncratic isolates shared a sensitivity pattern with lineage 58 (Fig. 4a, Extended Data Fig. 5), prompting a search for parallel mutations. We plotted the sensitivity phenotype against whole genome single nucleotide polymorphism trees constructed for each lineage. Each isolate on the tree was annotated with any vraFG and dlt operon mutations and/or genomic region gains or losses (colored by copy number relative to core genome). The sensitive phenotype (red/orange) correlated with vraF and vraG mutations in (A) lineage 20, (B) lineage 37, and (C) lineage 58. Presence or absence of genomic regions have no apparent trend in regards to the idiosyncratic sensitivity phenotype. While the vraFG mutants form a small fraction of their respective lineages, these isolates were found at multiple timepoints over six months apart, indicating that they are not a transient part of the population. (D) Since each represented a minority of their lineages, idiosyncratic isolates were excluded from other analyses presented in this study. However, on-person antagonisms remain depleted even if these isolates are retained (using simulations as in Fig. 3c in TSB). Significance P-value is derived from a two-sided permutation test with shuffled lineage labels (Methods), without multiple hypothesis correction. (E) Despite other sensitivities (Fig. 4d), the vraF-frameshift isolate (37.3) was less susceptible to phage infection than isolates with full-length vraF. However, we did not find this phage on subjects carrying lineage 37 (Methods). Data are presented as mean +/− standard deviation, n = 6, P values were calculated by unpaired t-test (two-sided). (F) Nevertheless, together with prior results from literature, these results are consistent with a conceptual model in which the vraFG complex could reduce cell envelope modifications, weakening resistance to cationic antimicrobials but providing protection against some phage. (G) Heatmap depicting the Minimum Inhibitory Concentration for antibiotics, lysis enzymes, and antimicrobial fatty acids, based on a screen with one replicate per condition. Polymyxin B (+5) and kanamycin (+4) are more positively charged than clindamycin (+1), ciprofloxacin (0), or cephalothin (−1). Idiosyncratic isolates 20.3 and 37.3 are more sensitive to cationic antimicrobials and lysostaphin than other isolates in their lineages. Both vraFG and dlt operon genotypes in lineage 58 share a similar MIC pattern to these idiosyncratic isolates, especially in response to Polymyxin B.

Figure 1:. High-throughput assay of interbacterial warfare from skin microbiomes of related people.

(A) Antagonism between different commensal strains of the same species (circles) could theoretically prevent transmission of certain strains despite frequent contact between individuals. Ecologically-relevant warfare would result in a depletion of antagonisms on individual people relative to random assortment. (B) We built a library of 122 Staphylococcus epidermidis isolates representing the strain diversity on healthy individuals from six families, and included 23 isolates from other skin species from these people. We screened every pairwise interaction (14,884 S. epidermidis intraspecies interactions plus 6,141 including other species) using high-throughput “spot-on-lawn” assays, in which producer cells are spotted on top of a lawn of receiver cells; we identified hundreds of cases of antagonism, which appear as Zones-Of-Inhibition (ZOIs; Methods). (C) Representative photograph of spot-on-lawn assay, using array layout 1. Several S. epidermidis spots produce antimicrobials that antagonize the S. epidermidis receiver lawn, resulting in ZOIs of a range of intensities. Some spots on the plate are blank or were excluded for poor growth (see Supplemental Table 1 for array information). (D) We developed a computer-aided image analysis pipeline to identify even subtle ZOIs that pass the inhibition threshold (red dash) and quantify the strength of antagonism, determined by the Area Under the Curve of the intensity traces (AUC, blue; Methods).

Figure 2:. Intra-species antagonism is common and does not segregate at higher phylogenetic levels.

(A) No evolutionary structure is observed in antimicrobial production (rows) or sensitivity (columns) beyond the lineage level, indicating that antimicrobial production and resistance are fast evolving and not conserved. The heatmap depicts the AUC for each ZOI calculated for pairwise interactions between skin isolates dereplicated to the lineage level (Methods). Each pair is screened twice, with each isolate serving as both spot and lawn. Inhibition of a receiver lawn by a producer spot is indicated by a lighter square, with intensity indicating the size of the ZOI. Lineages and species are sorted by phylogenetic similarity (see tree at left, not to scale beyond S. epidermidis), and dashed lines on the heatmap species boundaries. Extended Data Fig. 3 depicts this data without dereplicating S. epidermidis isolates into lineages. (B-C) The likelihood of observing antagonism between two lineages does not significantly differ between members of the same or different (B) species or (C) quorum sensing variant (agr type). Significance p-values are derived from two-sided permutation tests with shuffled group labels (Methods). (D) Conversely, intra-lineage antagonism is negligible compared to inter-lineage antagonism. (E) To assess the mechanisms of antagonistic interactions, we repeated the screen using cells and filtrates from 18 S. epidermidis producers against an array of bait lawns chosen based on interaction profiles (Extended Data Fig. 5). Filtrates were subjected to a variety of pre-treatments. Nine interactions that did not inhibit in the filtrate were retested on media containing excess iron to identify iron-suppressible interactions. The results suggest a large variety of interaction mechanisms which are not neatly predictable from genomic context in this screen (Extended Data Fig. 5).

Figure 3:. Antagonism is depleted among lineages co-residing on individuals.

We combined the antagonism screen results with isolate-inferred lineage compositions (Extended Data Fig. 2 and Supplementary Table 3) to study the ecological relevance of antagonism. (A) The abundance-weighted Antagonism Frequency (AF) of interactions between co-resident lineages (e.g. same person, same timepoint) is depleted relative to AF calculated across all subjects pooled together (red line). Each datapoint depicts the mean AF across samples from a subject. Codes indicate family and subject. (B) Lineage pairs in the same sample are less likely to antagonize each other than lineages from different samples. Significance P-values derived from two-sided permutation test (Methods) without multiple hypothesis correction. (C) Simulations of a null model that incorporates lineage composition on subjects as well as antagonism structure also support a significant depletion of antagonism in the observed AF (arrow) of coexisting lineages on individuals. In each simulation, the antagonism matrix is kept constant but lineages are shuffled across subjects. As each simulation has a different expected AF, we compare the mean difference between each subject’s AF (i.e. white datapoints from A) and the expected AF (i.e. red line from A) for observed and simulated data (mean ΔAF) Significance P-values derived from two-sided permutation test (Methods) without multiple hypothesis correction. (D) Filtered heatmap from Fig. 2a, depicting only lineages present on two selected families (see Supplementary Fig. 3 for all families). Intensity indicates AUC for the ZOI, and lineages are repeated if present on multiple subjects. (E) The Bray-Curtis Similarity between the lineage composition of sample pairs is negatively correlated with the AF calculated across the same samples (Spearman’s rank correlation, two-sided), supporting a role for antagonism in preventing transmission. One pair of samples was used for each subject pair. (F) Lineages that antagonize a higher proportion of other lineages reach higher abundance on individuals (Spearman’s rank correlation, two sided).

Figure 4:. On-person evolution drives sensitivity to multiple antimicrobials.

(A) Two lineages show significant intra-lineage variation, each containing isolates with sensitivity to the same producers, depicted using a section of the isolate-level heatmap (from Extended Data Fig. 3). Isolates from the same lineage are indicated by decimals and color. Intensity indicates AUC for the ZOI. The two idiosyncratic isolates depicted (denoted by *) have different sensitivity profiles from the rest of their lineages and also have loss-of-function mutations in the vraFG pathway, as labeled. (B) The phylogeny of isolates in lineage 37 indicates that the vraF frameshift and associated sensitivity was acquired, not ancestral. Lineage 20 has a similar pattern of recently acquired sensitivity (see Extended Data Fig. 10a). (C) By selecting on polymyxin B, we cultivated revertant isolate 37.3-r1 which corrected the vraF frameshift and has an ancestral resistance profile. Genome sequencing confirmed this reversion (Methods). (D) Isolates with vraFG mutations, but not their ancestors or the revertant, are sensitive to cationic antimicrobials, lysostaphin, and filtered supernatant from isolate 35.1. Minimum Inhibitory Concentration (MIC) assay results were normalized to isolate 35.1 and raw values are available in Extended Data Fig. 10g.