Propionibacterium acnes Strain Populations in the Human Skin Microbiome Associated with Acne

  • Author Footnotes
    8 These authors contributed equally to this work.
    Sorel Fitz-Gibbon
    Footnotes
    8 These authors contributed equally to this work.
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Author Footnotes
    8 These authors contributed equally to this work.
    Shuta Tomida
    Footnotes
    8 These authors contributed equally to this work.
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Bor-Han Chiu
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Lin Nguyen
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Christine Du
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA

    Department of Medicine, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Minghsun Liu
    Affiliations
    Department of Microbiology, Immunology and Molecular Genetics, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • David Elashoff
    Affiliations
    Department of Medicine, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Marie C. Erfe
    Affiliations
    Center for Immunotherapeutics Research, Los Angeles Biomedical Research Institute at Harbor-UCLA Medical Center, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Anya Loncaric
    Affiliations
    Department of Medicine, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Jenny Kim
    Affiliations
    Department of Medicine, David Geffen School of Medicine, UCLA, Los Angeles, California, USA

    Department of Dermatology, Veterans Affairs Greater Los Angeles Healthcare System, Los Angeles, California, USA
    Search for articles by this author
  • Robert L. Modlin
    Affiliations
    Department of Medicine, David Geffen School of Medicine, UCLA, Los Angeles, California, USA

    Department of Microbiology, Immunology and Molecular Genetics, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Jeff F. Miller
    Affiliations
    Department of Microbiology, Immunology and Molecular Genetics, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • Erica Sodergren
    Affiliations
    The Genome Institute at Washington University, St Louis, Missouri, USA
    Search for articles by this author
  • Noah Craft
    Affiliations
    Center for Immunotherapeutics Research, Los Angeles Biomedical Research Institute at Harbor-UCLA Medical Center, David Geffen School of Medicine, UCLA, Los Angeles, California, USA
    Search for articles by this author
  • George M. Weinstock
    Affiliations
    The Genome Institute at Washington University, St Louis, Missouri, USA
    Search for articles by this author
  • Huiying Li
    Correspondence
    4339 CNSI, 570 Westwood Plaza, Building 114, UCLA, Los Angeles, California 90095-1770, USA
    Affiliations
    Department of Molecular and Medical Pharmacology, Crump Institute for Molecular Imaging, David Geffen School of Medicine, UCLA, Los Angeles, California, USA

    UCLA-DOE Institute for Genomics and Proteomics, Los Angeles, California, USA
    Search for articles by this author
  • Author Footnotes
    8 These authors contributed equally to this work.
      The human skin microbiome has important roles in skin health and disease. However, bacterial population structure and diversity at the strain level is poorly understood. We compared the skin microbiome at the strain level and genome level of Propionibacterium acnes, a dominant skin commensal, between 49 acne patients and 52 healthy individuals by sampling the pilosebaceous units on their noses. Metagenomic analysis demonstrated that although the relative abundances of P. acnes were similar, the strain population structures were significantly different in the two cohorts. Certain strains were highly associated with acne, and other strains were enriched in healthy skin. By sequencing 66 previously unreported P. acnes strains and comparing 71 P. acnes genomes, we identified potential genetic determinants of various P. acnes strains in association with acne or health. Our analysis suggests that acquired DNA sequences and bacterial immune elements may have roles in determining virulence properties of P. acnes strains, and some could be future targets for therapeutic interventions. This study demonstrates a previously unreported paradigm of commensal strain populations that could explain the pathogenesis of human diseases. It underscores the importance of strain-level analysis of the human microbiome to define the role of commensals in health and disease.

      Abbreviations

      CRISPR
      Clustered Regularly Interspaced Short Palindromic Repeat
      HMP
      Human Microbiome Project
      rDNA
      ribosomal DNA
      RT
      ribotype

      Introduction

      The diversity of the human microbiota at the strain level and its association with human health and disease are largely unknown. However, many studies have shown that microbe-related human diseases are often caused by certain strains of a species, rather than the entire species being pathogenic. Examples include methicillin-resistant Staphylococcus aureus (
      • Chambers H.F.
      • Deleo F.R.
      Waves of resistance: Staphylococcus aureus in the antibiotic era.
      ;
      • Chen C.J.
      • Su L.H.
      • Lin T.Y.
      • et al.
      Molecular analysis of repeated methicillin-resistant Staphylococcus aureus infections in children.
      ;
      • Hansra N.K.
      • Shinkai K.
      Cutaneous community-acquired and hospital-acquired methicillin-resistant Staphylococcus aureus.
      ) and Escherichia coli O157 (
      • Tarr P.I.
      • Gordon C.A.
      • Chandler W.L.
      Shiga-toxin-producing Escherichia coli and haemolytic uraemic syndrome.
      ;
      • Chase-Topping M.
      • Gally D.
      • Low C.
      • et al.
      Super-shedding and the link between human infection and livestock carriage of Escherichia coli O157.
      ). Acne vulgaris (commonly called acne) is one of the most common skin diseases with a prevalence in up to 85% of teenagers and 11% of adults (
      • White G.M.
      Recent findings in the epidemiologic evidence, classification, and subtypes of acne vulgaris.
      ). Although the etiology and pathogenesis of acne are still unclear, microbial involvement is considered to be one of the main mechanisms contributing to the development of acne (
      • Cunliffe W.J.
      Looking back to the future - acne.
      ;
      • Bojar R.A.
      • Holland K.T.
      Acne and Propionibacterium acnes.
      ). In particular, Propionibacterium acnes has been hypothesized to be an important pathogenic factor (
      • Webster G.F.
      Inflammation in acne vulgaris.
      ). Antibiotic therapy targeting P. acnes has been a mainstay treatment for more than 30 years (
      • Leyden J.J.
      The evolving role of Propionibacterium acnes in acne.
      ). However, despite decades of study, it is still not clear how P. acnes contributes to acne pathogenesis while being a major commensal of the normal skin flora (
      • Gao Z.
      • Tseng C.H.
      • Pei Z.
      • et al.
      Molecular analysis of human forearm superficial skin bacterial biota.
      ;
      • Bek-Thomsen M.
      • Lomholt H.B.
      • Kilian M.
      Acne is not associated with yet-uncultured bacteria.
      ;
      • Cogen A.L.
      • Nizet V.
      • Gallo R.L.
      Skin microbiota: a source of disease or defence?.
      ;
      • Fierer N.
      • Hamady M.
      • Lauber C.L.
      • et al.
      The influence of sex, handedness, and washing on the diversity of hand surface bacteria.
      ;
      • Costello E.K.
      • Lauber C.L.
      • Hamady M.
      • et al.
      Bacterial community variation in human body habitats across space and time.
      ;
      • Grice E.A.
      • Kong H.H.
      • Conlan S.
      • et al.
      Topographical and temporal diversity of the human skin microbiome.
      ;
      • Dominguez-Bello M.G.
      • Costello E.K.
      • Contreras M.
      • et al.
      Delivery mode shapes the acquisition and structure of the initial microbiota across multiple body habitats in newborns.
      ). Whether P. acnes protects the human skin as a commensal bacterium or functions as a pathogenic factor in acne, or both, remains to be elucidated.
      Here we compared the skin microbiome at the strain level and genome level in 49 acne patients and 52 normal individuals using a combination of metagenomics and genome sequencing. First, for each sample, 16S ribosomal DNA (rDNA) was amplified, ∼400 clones were sequenced, and an average of 311 nearly full-length 16S rDNA sequences were analyzed. The population structure of P. acnes strains was determined in each sample. Second, each P. acnes strain was assigned an “acne index” by calculating its prevalence in acne patients based on the 16S rDNA metagenomic data. The P. acnes strains associated with the acne patient group were identified, as well as the strains enriched in the individuals with normal skin. This metagenomic approach is fundamentally different from prior approaches in determining disease associations; it is more powerful and less biased than traditional methods by bypassing the biases and selection in strain isolation and culturing. To our knowledge, this study has the largest number of individual skin microbiomes reported at the strain level to date. Finally, we sequenced 66 previously unreported P. acnes strains and compared 71 P. acnes genomes covering the major lineages of P. acnes found in the skin microbiota. By combining a metagenomic study of the skin microbiome and genome sequencing of this major skin commensal, this study provides insight into potential bacterial genetic determinants in acne pathogenesis and emphasizes the importance of strain-level analysis of the human microbiome to understand the role of commensals in health and disease.

      Results

       P. acnes dominates the pilosebaceous unit

      We characterized the microbiome in pilosebaceous units (“pores”) on the nose, collected from 49 acne patients and 52 individuals with normal skin. Nearly full-length 16S rDNA sequences were obtained using the Sanger method, which allowed us to analyze P. acnes at the strain level. After quality filtering, our final data set contained 31,461 16S rDNA sequences ranging from position 29 to position 1,483. Of the sequences, 27,358 matched to P. acnes with greater than 99% identity. Our data demonstrated that P. acnes dominates the microbiota of pilosebaceous units, accounting for 87% of the clones (Figure 1). Other commonly found species in pilosebaceous units included Staphylococcus epidermidis, Propionibacterium humerusii, and Propionibacterium granulosum, each representing 1–2.3% of the total clones. A total of 536 species-level operational taxonomic units belonging to 42 genera and six phyla were identified in the samples (Supplementary Table S1 online).
      Figure thumbnail gr1
      Figure 1Propionibacterium acnes was dominant in pilosebaceous units in both acne patients and individuals with normal skin. By 16S ribosomal DNA (rDNA) sequencing, P. acnes sequences accounted for 87% of all the clones. Species with a relative abundance of >0.35% are listed in order of relative abundance. Species distribution from a metagenomic shotgun sequencing of pooled samples from normal individuals confirmed the high abundance of P. acnes in pilosebaceous units, as shown in the far right column.
      To bypass the potential biases due to PCR amplification and due to uneven numbers of 16S rDNA gene copies among different species, we performed a metagenomic shotgun sequencing of the total DNA pooled from the pilosebaceous unit samples of 22 additional normal individuals. Microbial species were identified by mapping metagenomic sequences to reference genomes. The results confirmed that P. acnes was the most abundant species (89%) (Figure 1). This is consistent with the results obtained from 16S rDNA sequencing (87%).

       Different P. acnes strain populations in acne

      There was no statistically significant difference in the relative abundance of P. acnes when comparing acne patients and normal individuals. We next examined whether there were differences at the strain level of P. acnes by extensively analyzing the P. acnes 16S rDNA sequences. We define each unique 16S rDNA sequence as a 16S rDNA allele type, called a ribotype (RT). The most abundant P. acnes sequence was defined as ribotype 1 (RT1); all other defined ribotypes have ≥99% sequence identity to RT1. Similar to the distributions seen at higher taxonomical levels (
      • Bik E.M.
      • Long C.D.
      • Armitage G.C.
      • et al.
      Bacterial diversity in the oral cavity of 10 healthy individuals.
      ), at the strain level a few ribotypes were highly abundant in the samples with a significant number of rare ribotypes (Supplementary Figure S1 online). After careful examination of the sequence chromatograms and manual correction of the sequences, a total of 11,009 ribotypes were assigned to the P. acnes 16S rDNA sequences. Most of the minor ribotypes were singletons. On average, each individual harbored 3±2 P. acnes ribotypes with three or more clones. On the basis of the genome sequences described below, all the sequenced P. acnes strains have three identical copies of 16S rDNA genes (note in Supplementary Information online). This allowed us to compare the P. acnes strain populations in individuals based on the 16S rDNA sequences. The top 10 major ribotypes with more than 60 clones and found in multiple subjects are shown in Table 1.
      Table 1Top 10 most abundant ribotypes found in pilosebaceous units
      Ribotype (RT)Nucleotide changes from RT1Number of subjectsNumber of clonesPercentage of clones from acne patients
      1 The percentage was calculated after the number of clones of each ribotype was normalized by the total number of clones in acne patients (acne index).
      Percentage of clones from normal individuals
      2 The percentage was calculated after the number of clones of each ribotype was normalized by the total number of clones in normal individuals.
      P-value
      3 Mann–Whitney–Wilcoxon rank-sum test.
      RT1905,53648%52%0.84
      RT2T854C481,21351%49%0.36
      RT3T1007C602,10440%60%0.092
      RT4G1058C, A1201C2327584%16%0.049
      RT5G1058C1520599%1%0.00050
      RT6T854C, C1336T112621%99%0.025
      RT7G529A1018899%1%0.12
      RT8G1004A, T1007C5239100%0%0.024
      RT9G1268A46899%1%0.29
      RT10T554C, G1058C561100%0%0.024
      1 The percentage was calculated after the number of clones of each ribotype was normalized by the total number of clones in acne patients (acne index).
      2 The percentage was calculated after the number of clones of each ribotype was normalized by the total number of clones in normal individuals.
      3 Mann–Whitney–Wilcoxon rank-sum test.
      Analysis of the top 10 ribotypes showed both disease-specific and health-specific associations. The three most abundant ribotypes (RT1, RT2, and RT3) were fairly evenly distributed among acne and normal individuals. However, the next seven major ribotypes were significantly skewed in their distributions (Table 1). Ribotypes 4, 5, 7, 8, 9, and 10 were found predominantly in acne patients, with four of these six ribotypes being statistically significantly enriched in acne (P<0.05, Wilcoxon test). Ribotypes 4, 5, and 10 contain a nucleotide substitution G1058C in the 16S rDNA sequences, which has previously been shown to confer increased resistance to tetracycline (
      • Ross J.I.
      • Eady E.A.
      • Cove J.H.
      • et al.
      16S rRNA mutation associated with tetracycline resistance in a gram-positive bacterium.
      ,
      • Ross J.I.
      • Snelling A.M.
      • Eady E.A.
      • et al.
      Phenotypic and genotypic characterization of antibiotic-resistant Propionibacterium acnes isolated from acne patients attending dermatology clinics in Europe, the U.S.A., Japan and Australia.
      ). However, only a small percentage of the subjects in our study harboring these ribotypes had been treated with antibiotics (Supplementary Table S2 online); therefore, enrichment of these three ribotypes in the acne group was not correlated with antibiotic treatment. This is consistent with previous studies, which showed that previous use of antibiotics was not always associated with the presence of antibiotic-resistant strains and that some patients who were not previously treated with antibiotics harbored strains already resistant to antibiotics (
      • Dreno B.
      • Reynaud A.
      • Moyse D.
      • et al.
      Erythromycin-resistance of cutaneous bacterial flora in acne.
      ;
      • Coates P.
      • Vyakrnam S.
      • Eady E.A.
      • et al.
      Prevalence of antibiotic-resistant propionibacteria on the skin of acne patients: 10-year surveillance data and snapshot distribution study.
      ). On the other hand, one ribotype, RT6, although detected in only 11 subjects, was strongly associated with normal skin (P=0.025, Wilcoxon test; Table 1). Its relative abundance in our normal group was similar to that found in the healthy cohort data from the Human Microbiome Project (HMP) (see Supplementary Information and Supplementary Figure S2 online). The percentage of positive subjects (11/52) was similar as well. Three of the 14 HMP subjects had RT6 found in the anterior nares, and one additional subject had RT6 in the left retroauricular crease.
      On the basis of the distributions of the top 10 ribotypes, statistical analysis using several different tests showed significant differences in P. acnes population structure between acne and normal skin (Supplementary Table S3 online). This is consistent with a principal coordinate analysis, in which acne samples and normal skin samples were separated mostly by principal coordinates 1 and 2 (Supplementary Figure S3 online), explaining 44% and 20% of the variation, respectively.
      To examine whether different individuals share similar P. acnes population structures, we clustered the samples on the basis of the relative abundances of the top 10 ribotypes. Five main microbiome types were observed at the P. acnes strain level (microbiome types I to V). Types IV and V, which are dominated by P. acnes RT4 and RT5, respectively, were mainly found in acne patients (Figure 2 and Supplementary Figure S4 online). The same five main microbiome types were observed in the HMP data and the data from
      • Grice E.A.
      • Kong H.H.
      • Conlan S.
      • et al.
      Topographical and temporal diversity of the human skin microbiome.
      (see Supplementary Information and Supplementary Figure S5 online).
      Figure thumbnail gr2
      Figure 2Distribution of the top 10 most abundant Propionibacterium acnes ribotypes in acne patients and individuals with normal skin. Each column represents the percentage of the top 10 ribotypes identified in each subject. The average P. acnes clone number per subject was 262, and the average clone number of the top 10 ribotypes was 100. Five major microbiome types at the P. acnes strain level were observed in the data. Types IV and V were mostly found in acne patients. Two samples (one from acne and one from normal skin) with <50 P. acnes 16S ribosomal DNA (rDNA) sequences are not displayed.

       Genome sequence analysis of 71 P. acnes strains

      All of the top 10 most abundant ribotypes differ from RT1 by only one or two nucleotide changes in the 16S rDNA sequence (Table 1). To determine whether such small changes in the 16S rDNA sequence reflect the lineages and evolutionary history at the genome level, we selected 66 P. acnes isolates representing major ribotypes 1, 2, 3, 4, 5, 6, and 8, as well as two minor ribotypes, 16 and 532, for genome sequencing. The genomes of these 66 isolates were fully sequenced and assembled to high-quality drafts or complete genomes with ≥50X coverage. Five other P. acnes genomes, KPA171202 (
      • Bruggemann H.
      • Henne A.
      • Hoster F.
      • et al.
      The complete genome sequence of Propionibacterium acnes, a commensal of human skin.
      ), J165, J139, SK137, and SK187, were publicly available and were included in our analysis. We constructed a phylogenetic tree based on 96,887 unique single-nucleotide polymorphism positions in the core genome obtained from these 71 P. acnes genomes. Most of the genomes with the same ribotypes clustered together. The tree suggests that the 16S rDNA ribotypes do represent the relationship of the lineages to a large extent, and that16S rDNA sequence is a useful molecular marker to distinguish major P. acnes lineages (Figure 3 and Supplementary Figure S6 online).
      Figure thumbnail gr3
      Figure 3Genome comparison of 71 Propionibacterium acnes strains showed that the genomes of ribotypes 4 and 5 (RT4 and RT5) are distinct from those of others. Two chromosomal regions, loci 1 and 2, are unique to clade IA-2 and one other genome HL086PA1. Clade IA-2 consists of mainly RT4 and RT5 that were highly enriched in acne. The presence of a plasmid (locus 3) is also characteristic of RT4 and RT5. Each row represents a P. acnes genome colored according to the ribotypes. Rows are ordered by the phylogeny calculated based on the single-nucleotide polymorphisms (SNPs) in the P. acnes core genome. Only the topology is shown. The clades were named based on their recA types (IA, IB, and II). Columns represent predicted open reading frames (ORFs) in the genomes and are ordered by ORF positions along the finished genome HL096PA1, which encodes a 55-Kb plasmid. Only the first 300 ORFs on the chromosome (on the left) and all the ORFs on the plasmid (on the right) are shown. The colored plasmid regions represent genes on contigs that match exclusively to the HL096PA1 plasmid region. The genes that fall on contigs that clearly extend beyond the plasmid region are likely to be chromosomally located and are colored in gray. Acne index for the ribotypes was calculated based on the percentage of clones of each ribotype found in acne as shown in column 5 in .

       Genetic elements detected in P. acnes

      We further performed comparative genome analysis among all 71 genomes grouped by ribotypes. Our analysis revealed potential genetic elements by which acne-associated strains could contribute to acne pathogenesis and the elements by which health-associated strains could contribute to maintaining skin health. Specifically, we describe here the unique genome regions of RT4 and RT5, which had a strong association with acne, and RT6, which was found to be enriched in normal skin. Three distinct regions, loci 1, 2, and 3, were found almost exclusively in strains that belong to clade IA-2 in the phylogenetic tree. Clade IA-2 consists of mainly RT4 and RT5 (Figure 3 and Supplementary Figure S7 online). Loci 1 and 2 are located on the chromosome. Locus 1 contains prophage-related genes and appears to be a genomic island. Locus 2 has plasmid integration sites and thus could be derived from a plasmid sequence. Locus 3 appears to be on a large mobile genetic element, likely a plasmid. The plasmid is ∼55Kb long and has inverted terminal repeats according to our finished genome HL096PA1 (Supplementary Information online). The sequence data suggest that the plasmid is linear and possibly originated from a phage (
      • Hinnebusch J.
      • Tilly K.
      Linear plasmids and chromosomes in bacteria.
      ). All but 1 of the 15 genomes of RT4 and RT5 have at least 60% of the genes of the plasmid represented, and all of them have regions homologous to the inverted terminal repeats in the plasmid, suggesting that they harbor the same or a similar linear plasmid (Figure 3). The copy number of the plasmid in the genomes ranges from 1 to 3 based on genome sequencing coverage, which was confirmed by quantitative PCR (Supplementary Figures S8 and S9 online).
      The fact that acne-enriched RT4 and RT5 strains carry a linear plasmid and two unique loci of genomic islands suggests that these plasmid and chromosomal regions may have a role in acne pathogenesis. In fact, the linear plasmid encodes a tight adhesion (Tad) locus, which has been suggested to have a role in virulence in other organisms (
      • Kachlany S.C.
      • Planet P.J.
      • Bhattacharjee M.K.
      • et al.
      Nonspecific adherence by Actinobacillus actinomycetemcomitans requires genes widespread in bacteria and archaea.
      ;
      • Schreiner H.C.
      • Sinatra K.
      • Kaplan J.B.
      • et al.
      Tight-adherence genes of Actinobacillus actinomycetemcomitans are required for virulence in a rat model.
      ). The complete Tad locus is found in all but 1 of the 15 genomes of RT4 and RT5, and is only occasionally found in other ribotypes. In addition, in locus 2, a Sag gene cluster is encoded, which has been shown to contribute to hemolytic activity in pathogens (
      • Nizet V.
      • Beall B.
      • Bast D.J.
      • et al.
      Genetic locus for streptolysin S production by group A streptococcus.
      ;
      • Fuller J.D.
      • Camus A.C.
      • Duncan C.L.
      • et al.
      Identification of a streptolysin S-associated gene cluster and its role in the pathogenesis of Streptococcus iniae disease.
      ;
      • Humar D.
      • Datta V.
      • Bast D.J.
      • et al.
      Streptolysin S and necrotising infections produced by group G streptococcus.
      ). Supplementary Table S4 online summarizes the genes that are mostly unique to RT4 and RT5, several of which have essential roles in virulence in other organisms. We speculate that some of these genes encoded in RT4 and RT5 may increase virulence, promote stronger adherence to the human host, or induce a pathogenic host immune response.
      In genome comparison analysis, we found that all the genomes of RT2 and RT6 encode CRISPRs (Clustered Regularly Interspaced Short Palindromic Repeats). Among the sequenced genomes, RT2 and RT6 are the only ribotypes encoding CRISPRs. CRISPRs have been shown to confer protective “immunity” against viruses, phages, and plasmids (
      • Horvath P.
      • Barrangou R.
      CRISPR/Cas, the immune system of bacteria and archaea.
      ;
      • Makarova K.S.
      • Haft D.H.
      • Barrangou R.
      • et al.
      Evolution and classification of the CRISPR-Cas systems.
      ). The CRISPR locus encoded in P. acnes consists of a series of cas genes—cas3, cse1, cse2, cse4, cas5e, cse3, cas1, and cas2—which are homologous to the CRISPR locus reported in E. coli (Supplementary Figure S10 online) and the CRISPR4 locus in Streptococcus thermophilus (
      • Horvath P.
      • Barrangou R.
      CRISPR/Cas, the immune system of bacteria and archaea.
      ).
      CRISPR arrays are composed of a cluster of identical repetitive sequences separated by spacer sequences of similar length but with different nucleotide sequences. It has been found that spacer sequences are identical, or with one or two mismatches, to phage or plasmid DNA sequences. A total of 39 spacer sequences were found in 8 P. acnes strains, 25 of which were unique, as shown in Table 2. As expected, most of the identifiable spacers target to known P. acnes phage sequences. However, among the unique CRISPR spacer sequences, one matched locus 2 on the chromosome and three matched the plasmid region (locus 3) in P. acnes genomes of mainly RT4 and RT5. This suggests that these loci may have been acquired by RT4 and RT5 strains, whereas the genomes of RT2 and RT6 may be capable of protecting against the invasion of the plasmids or other foreign DNA through the CRISPR mechanism.
      Table 2CRISPR spacer sequences found in the genomes of RT2 and RT6
      RibotypeStrainSpacer numberSpacer sequenceBLAST resultMatch found
      RT2HL001PA11CATGGCCTGCACACCAGGCGCTTTTAGCACCTNo hits
      2CATGGCCTGCACACCAGGCGCTTTTAGCACCTNo hits
      3CATGGCCTGCACACCAGGCGCTTTTAGCACCTNo hits
      4GGCGTATGACGAGTTGTGGTCGGCGTTTCCTCP. acnes phage PA6 gp15 (minor tail protein)
      5CGGTGTTAACGGCTTGCCTGGCTTGGATGGAGNo hits
      RT2HL060PA11CGCCTACCGTCAGCTGACTCACGCCTCCGCGTTNo hits
      2TCACACCAGTCATCAGCGTCATAGTCCTCTCGGNo hits
      RT2HL082PA21GGCTCAGCCCTGCCCGATGCCTACGCCAAATGGC. leptum DSM 753 CLOLEP_00129 (cell wall–associated hydrolases (invasion-associated proteins))Locus 3
      2TCACACCAGTCATCAGCGTCATAGTCCTCTCGGNo hits
      RT2HL103PA11CACCGGGCCCATCCCGGTCGGCCTCCTGAAAGGC. leptum DSM 753 CLOLEP_00135Locus 3
      RT2HL106PA11GATCGAGTTGGCTGAGTCGAAGGTGTTGCGGTTP. acnes phage PA6 gp16 (conserved protein)

      P. acnes phage PAD20 gp16
      2CTGCTCATCGCTCAGCTCCTGCGCCTCATCACANo hits
      3CTGCGCCAACAGCCGCATCTGATCCGAATACGGP. acnes phage PA6 gp3 (phage portal protein)
      4CGCAGCAATCTCAGAAGGCCACAACAAGTTCGTP. acnes phage PA6 gp7 (conserved protein)

      P. acnes phage PAD20 gp7

      P. acnes phage PAS50 gp7
      5CAAATCACCCAAGCCCAACACGCCGCCACCACCNo hits
      6TGTCACCGATTCAATGTATCTATGAGTGGTGTANo hits
      7TTGGGTGGGTGAGGTCGGGTCGTCAGTCATGAGNo hits
      8GTCGATGTCGAGATTGGCCTGGGGGTCCATGTCC. leptum DSM 753 CLOLEP_00142Locus 3
      9ACGTCGTGAACGTACCCCTTGACGGAGACGGCANo hits
      RT2J1391CGAGGGCTACCACGTGGTCGATTTGGACTGTCGC. leptum DSM 753 CLOLEP_00167

      P. acnes SK137 HMPREF0675_3193 (domain of unknown function)
      Locus 2
      2CAGGCGCTCCACTCCCTCGCCCTGGCCACCAACNo hits
      RT6HL110PA3HL110PA41CTATGTGGACAGTGTTGGTTACTGTGGGGGGAAP. acnes phage PA6 intergenic region between gp45 and gp46
      2GCACTGCACCGATATCGTCGTGGCTGTCACTTGNo hits
      3CCCAGACAACCTCGACAACCTGTTCAGGGGATGP. acnes phage PAS50 gp25
      4CATGGCTAGCCCGGATTTTTGGCTGCCTGAGCGP. acnes phage PA6 gp34 (multidrug resistance protein-like transporters)

      P. acnes phage PAD20 gp34 (DNA helicase)
      5CGGCCTGCGGCAGATTTTTGTTGCGTTGAATCCP. acnes phage PA6 gp14 (tape measure protein)

      P. acnes phage PAD20 gp14 (tape measure protein)

      P. acnes phage PAS50 gp14 (tape measure protein)
      6CGGGCAGAGGATGTGTTGCTCGTTCCTGGATGGP. acnes phage PA6 gp32 (CHC2 zinc-finger)

      P. acnes phage PAD20 gp32 (DNA primase)

      P. acnes phage PAS50 gp32 (DNA primase)
      7GTTACGCTGGAACCCCCAATGAACACGCGAGAAP. acnes phage PAD42 major head protein gene

      P. acnes phage PAD20 major head protein gene

      P. acnes phage PAD9 major head protein gene

      P. acnes phage PAS40 major head protein gene

      P. acnes phage PAS12 major head protein gene

      P. acnes phage PAS10 major head protein gene

      P. acnes phage PAD21 major head protein gene

      P. acnes phage PAS2 major head protein gene

      P. acnes phage PA6 gp6 (phage capsid family)

      P. acnes phage PAS50 gp6 major head protein gene
      8CGAGGGCTACCACGTGGTCGATTTGGACTGTCGC. leptum DSM 753 CLOLEP_00167

      P. acnes SK137 HMPREF0675_3193 (domain of unknown function)
      Locus 2
      9CAGGCGCTCCACTCCCTCGCCCTGGCCACCAACNo hits
      Abbreviations: BLAST, Basic Local Alignment Search Tool; CRISPR, Clustered Regularly Interspaced Short Palindromic Repeat; C. leptum, Clostridium leptum; P. acnes, Propionibacterium acnes; RT, ribotype.

      Discussion

      Our study of the human skin microbiome associated with acne provides a previously unreported portrait of the microbiota of pilosebaceous units at the bacterial strain level. As P. acnes is the major skin commensal bacterium found in both acne and healthy skin, this strain-level analysis is important to help understand the role of P. acnes in acne pathogenesis and in skin health. We demonstrate a strong association between strains of RT4 and RT5 with acne and a strong association between strains of RT6 and healthy skin, each with unique genetic elements. Other P. acnes strains, including ribotypes 7, 8, 9, and 10, or interactions among different strains, may also contribute to the development of the disease. In addition, host factors, such as hormone level, sebum production, and physical changes in the pilosebaceous unit, may also have a role in acne pathogenesis. Further studies aimed at identifying the specific functions of these strains, host factors in the development of acne, as well as the associations of microbiome characteristics with the subtypes of acne (comedonal, pustular, inflammatory, cystic, etc.) with larger cohort sizes, may improve our understanding of the molecular mechanisms of the disease. These studies may also help develop a targeted therapeutic approach to treat this extremely common and sometimes disfiguring skin disease.
      Our metagenomic approach in revealing the association of P. acnes strains with the disease or health is, to our knowledge, previously unreported, and is more powerful than previous studies using traditional methods (
      • Lomholt H.B.
      • Kilian M.
      Population genetic analysis of Propionibacterium acnes identifies a subpopulation and epidemic clones associated with acne.
      ;
      • McDowell A.
      • Gao A.
      • Barnard E.
      • et al.
      A novel multilocus sequence typing scheme for the opportunistic pathogen Propionibacterium acnes and characterization of type I cell surface-associated antigens.
      ). Because the skin microbiota of each individual and each skin site may harbor “good,” “neutral,” and “bad” strains at the same time, which may have different growth rates under in vitro culturing conditions, culturing a few isolates from a disease lesion or healthy skin site may not provide an accurate and unbiased measurement of the association of the strains with the disease or health. The sampling technique and disease associations in this study do not depend on sampling locations, on the presence of lesions in the sampling field, or on inherently biased culture techniques. Although sampling lesional skin intentionally may yield interesting results, these results would not be capable of defining the disease associations that unbiased sampling can. The metagenomic approach used in this study to identify underlying strain differences in acne might also be applied to the study of other disease/health associations with commensal or pathogenic bacteria. Ultimately, these studies could lead to targeted therapeutics to restore natural commensal population structures, and could help determine whether therapeutic modulations of the microbiota can return the host to a state of health.

      Materials and Methods

       Subjects

      Subjects with acne and subjects with normal skin were recruited from various clinics in Southern California, including private practice, managed care, and public hospital settings, as well as outside of dermatology clinics, to best represent the diversity of populations and history of medical care. The subject data are available at dbGaP (http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000263.v1.p1). The diagnosis of acne was made by board-certified dermatologists. The presence of acne was graded on a scale of 0 to 5 relating closely to the Global Acne Severity Scale (
      • Dreno B.
      • Poli F.
      • Pawin H.
      • et al.
      Development and evaluation of a Global Acne Severity Scale (GEA Scale) suitable for France and Europe.
      ). Grades were recorded for both the face and the nose separately, where 0 represents normal skin and 5 represents the most severe inflammatory cystic acne. In acne patients, the grades of the face ranged from 1 to 5, with an average of 2.1, and the grades of the nose ranged from 0 to 2, with an average of 0.3. The presence of scarring was also noted. Subjects with normal skin were determined by board-certified dermatologists and were defined as people who had no acneiform lesions on the face, chest, or back. They were also excluded if they had other skin problems that the investigators felt would affect sampling or the microbial population on the skin. Among the 101 subjects, 59 were female (31 acne patients and 28 normal subjects) and 42 were male (18 acne patients and 24 normal subjects). The average age of the acne cohort was 22.2, and the average age of the normal cohort was 29.6. There was no significant difference in ethnicity between the acne and normal populations. The subjects responded to a written questionnaire, administered by a physician or a well-trained study coordinator who went over each question with the subjects. Most of the subjects had not been treated for acne in the past or were not being treated when samples were collected (Supplementary Table S2). Only 9 out of 78 subjects, who provided treatment information, were being treated for acne when samples were taken. Among the nine subjects, two were being treated with antibiotics, five were being treated with topical retinoids, one was being treated with both antibiotics and retinoids, and one did not list the treatment. We also asked subjects for acne treatment history in the past (anytime in their life). Out of 73 subjects who provided treatment history, 18 had been treated for acne in the past. Among them, seven had been treated with antibiotics, eight had been treated with retinoids, two had been treated with both antibiotics and retinoids, and one did not list the treatment. All subjects provided written informed consent. All protocols and consent forms were approved by the institutional review boards of both the University of California, Los Angeles and the Los Angeles Biomedical Research Institute. The study was conducted in adherence to the Declaration of Helsinki Principles.

       Samples

      Skin microcomedone (white head or black head) samples were taken from the nose of the subjects using Bioré Deep Cleansing Pore Strips (Kao Brands Company, Cincinnati, OH) by following the instructions of the manufacturer. Clean gloves were used for each sampling. After being removed from the nose, the strip was immediately placed into a 50-ml sterile tube and kept on ice or at 4°C. The cells were lysed within 4hours in most of the cases.

       Metagenomic DNA extraction, 16S rDNA amplification, cloning, and sequencing

      Individual microcomedones were isolated from the adhesive nose strip using sterile forceps. Genomic DNA was extracted using the QIAamp DNA Micro Kit (Qiagen, Valencia, CA). 16S rDNA was amplified and cloned according to the protocol by HMP, which is described in detail in Supplementary Information online. Nearly full-length sequences were obtained by the Sanger method.

       16S rDNA sequence analysis

      Base calling and quality was determined with Phred (
      • Ewing B.
      • Green P.
      Base-calling of automated sequencer traces using phred. II. Error probabilities.
      ;
      • Ewing B.
      • Hillier L.
      • Wendl M.C.
      • et al.
      Base-calling of automated sequencer traces using phred. I. Accuracy assessment.
      ). Bidirectional reads were assembled and aligned to a core set of NAST-formatted sequences (rRNA16S.gold) using AmosCmp16Spipeline and NAST-ier (
      • Haas B.J.
      • Gevers D.
      • Earl A.M.
      • et al.
      Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.
      ). Suspected chimeras were identified using ChimeraSlayer and WigeoN (
      • Haas B.J.
      • Gevers D.
      • Earl A.M.
      • et al.
      Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.
      ). 16S rDNA sequences were extensively manually examined. Chromatograms were visually inspected at all bases, with a Phred quality score of <30. Appropriate corrections were applied. QIIME (
      • Caporaso J.G.
      • Kuczynski J.
      • Stombaugh J.
      • et al.
      QIIME allows analysis of high-throughput community sequencing data.
      ) was used to cluster the sequences into operational taxonomic units.

       P. acnes isolation and genotyping

      Colonies with the macroscopic characteristics of P. acnes were picked from each sample plate and were passed twice. The ribotype of each isolate was determined by PCR amplification and sequencing of the full length of the 16S rDNA gene by the Sanger method.

       Whole-genome shotgun sequencing, assembly, and annotation

      Genome HL096PA1 was sequenced using Roche/454 FLX and was assembled using a combination of PHRAP/CONSED (
      • Gordon D.
      • Abajian C.
      • Green P.
      Consed: a graphical tool for sequence finishing.
      ) and GSMAPPER (Roche, Branford, CT) with extensive manual editing in CONSED. The remaining 65 genomes were sequenced using Illumina/Solexa GAIIx (Illumina, San Diego, CA). Sequence data sets were processed by quality trimming and were assembled using Velvet (
      • Zerbino D.R.
      • Birney E.
      Velvet: algorithms for de novo short read assembly using de Bruijn graphs.
      ). Coding sequences were predicted using GeneMark (
      • Borodovsky M.
      • McIninch J.
      Recognition of genes in DNA sequence with ambiguities.
      ) and GLIMMER. The final gene set was processed through a suite of protein categorization tools consisting of Interpro, psort-b, and KEGG. A more detailed protocol can be found at http://hmpdacc.org/doc/sops/reference_genomes/annotation/WUGC_SOP_DACC.pdf.

       Comparative genome analysis

      A total of 71 P. acnes genome sequences were compared using Nucmer (
      • Kurtz S.
      • Phillippy A.
      • Delcher A.L.
      • et al.
      Versatile and open software for comparing large genomes.
      ). Phylogenetic analysis was performed using MEGA5 (
      • Tamura K.
      • Dudley J.
      • Nei M.
      • et al.
      MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0.
      ). CRISPRFinder (
      • Grissa I.
      • Vergnaud G.
      • Pourcel C.
      CRISPRFinder: a web tool to identify clustered regularly interspaced short palindromic repeats.
      ) was used to identify the CRISPR repeat-spacer sequences.
      The data reported in this paper are listed in Supplementary Information online and archived at GenBank.

      ACKNOWLEDGMENTS

      We thank G Kasimatis, B Shi, EE Curd, R Yan, M Wong, and J Liu for comments and technical support. We thank C Lee for performing statistical analyses in the initial phase. We also thank Z Guo and CS Miller for critical reading of the manuscript. This research was funded as one of the Demonstration Projects by the NIH Human Microbiome Project (HMP). It was supported by grant UH2AR057503 from NIH/NIAMS and grant U54HG004968 from NIH.
      Author contributions
      SF-G and ST analyzed the data; B-HC, LN, and CD performed experiments; DE performed some of the initial statistical analyses; MCE, AL, JK, RLM, and NC collected samples; ES and GMW directed sequencing, genome assembly, and annotation; HL, ML, RLM, and JFM conceived the demonstration project in the initial phase; HL designed and directed the project, analyzed the data, and wrote the paper; and ST, SF-G, and NC co-wrote the paper.

      SUPPLEMENTARY MATERIAL

      Supplementary material is linked to the online version of the paper at http://www.nature.com/jid

      REFERENCES

        • Bek-Thomsen M.
        • Lomholt H.B.
        • Kilian M.
        Acne is not associated with yet-uncultured bacteria.
        J Clin Microbiol. 2008; 46: 3355-3360
        • Bik E.M.
        • Long C.D.
        • Armitage G.C.
        • et al.
        Bacterial diversity in the oral cavity of 10 healthy individuals.
        ISME J. 2010; 4: 962-974
        • Bojar R.A.
        • Holland K.T.
        Acne and Propionibacterium acnes.
        Clin Dermatol. 2004; 22: 375-379
        • Borodovsky M.
        • McIninch J.
        Recognition of genes in DNA sequence with ambiguities.
        Biosystems. 1993; 30: 161-171
        • Bruggemann H.
        • Henne A.
        • Hoster F.
        • et al.
        The complete genome sequence of Propionibacterium acnes, a commensal of human skin.
        Science. 2004; 305: 671-673
        • Caporaso J.G.
        • Kuczynski J.
        • Stombaugh J.
        • et al.
        QIIME allows analysis of high-throughput community sequencing data.
        Nat Methods. 2010; 7: 335-336
        • Chambers H.F.
        • Deleo F.R.
        Waves of resistance: Staphylococcus aureus in the antibiotic era.
        Nat Rev Microbiol. 2009; 7: 629-641
        • Chase-Topping M.
        • Gally D.
        • Low C.
        • et al.
        Super-shedding and the link between human infection and livestock carriage of Escherichia coli O157.
        Nat Rev Microbiol. 2008; 6: 904-912
        • Chen C.J.
        • Su L.H.
        • Lin T.Y.
        • et al.
        Molecular analysis of repeated methicillin-resistant Staphylococcus aureus infections in children.
        PLoS One. 2010; 5: e14431
        • Coates P.
        • Vyakrnam S.
        • Eady E.A.
        • et al.
        Prevalence of antibiotic-resistant propionibacteria on the skin of acne patients: 10-year surveillance data and snapshot distribution study.
        Br J Dermatol. 2002; 146: 840-848
        • Cogen A.L.
        • Nizet V.
        • Gallo R.L.
        Skin microbiota: a source of disease or defence?.
        Br J Dermatol. 2008; 158: 442-455
        • Costello E.K.
        • Lauber C.L.
        • Hamady M.
        • et al.
        Bacterial community variation in human body habitats across space and time.
        Science. 2009; 326: 1694-1697
        • Cunliffe W.J.
        Looking back to the future - acne.
        Dermatology. 2002; 204: 167-172
        • Dominguez-Bello M.G.
        • Costello E.K.
        • Contreras M.
        • et al.
        Delivery mode shapes the acquisition and structure of the initial microbiota across multiple body habitats in newborns.
        Proc Natl Acad Sci USA. 2010; 107: 11971-11975
        • Dreno B.
        • Poli F.
        • Pawin H.
        • et al.
        Development and evaluation of a Global Acne Severity Scale (GEA Scale) suitable for France and Europe.
        J Eur Acad Dermatol Venereol. 2011; 25: 43-48
        • Dreno B.
        • Reynaud A.
        • Moyse D.
        • et al.
        Erythromycin-resistance of cutaneous bacterial flora in acne.
        Eur J Dermatol. 2001; 11: 549-553
        • Ewing B.
        • Green P.
        Base-calling of automated sequencer traces using phred. II. Error probabilities.
        Genome Res. 1998; 8: 186-194
        • Ewing B.
        • Hillier L.
        • Wendl M.C.
        • et al.
        Base-calling of automated sequencer traces using phred. I. Accuracy assessment.
        Genome Res. 1998; 8: 175-185
        • Fierer N.
        • Hamady M.
        • Lauber C.L.
        • et al.
        The influence of sex, handedness, and washing on the diversity of hand surface bacteria.
        Proc Natl Acad Sci USA. 2008; 105: 17994-17999
        • Fuller J.D.
        • Camus A.C.
        • Duncan C.L.
        • et al.
        Identification of a streptolysin S-associated gene cluster and its role in the pathogenesis of Streptococcus iniae disease.
        Infect Immun. 2002; 70: 5730-5739
        • Gao Z.
        • Tseng C.H.
        • Pei Z.
        • et al.
        Molecular analysis of human forearm superficial skin bacterial biota.
        Proc Natl Acad Sci USA. 2007; 104: 2927-2932
        • Gordon D.
        • Abajian C.
        • Green P.
        Consed: a graphical tool for sequence finishing.
        Genome Res. 1998; 8: 195-202
        • Grice E.A.
        • Kong H.H.
        • Conlan S.
        • et al.
        Topographical and temporal diversity of the human skin microbiome.
        Science. 2009; 324: 1190-1192
        • Grissa I.
        • Vergnaud G.
        • Pourcel C.
        CRISPRFinder: a web tool to identify clustered regularly interspaced short palindromic repeats.
        Nucleic Acids Res. 2007; 35: W52-W57
        • Haas B.J.
        • Gevers D.
        • Earl A.M.
        • et al.
        Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.
        Genome Res. 2011; 21: 494-504
        • Hansra N.K.
        • Shinkai K.
        Cutaneous community-acquired and hospital-acquired methicillin-resistant Staphylococcus aureus.
        Dermatol Ther. 2011; 24: 263-272
        • Hinnebusch J.
        • Tilly K.
        Linear plasmids and chromosomes in bacteria.
        Mol Microbiol. 1993; 10: 917-922
        • Horvath P.
        • Barrangou R.
        CRISPR/Cas, the immune system of bacteria and archaea.
        Science. 2010; 327: 167-170
        • Humar D.
        • Datta V.
        • Bast D.J.
        • et al.
        Streptolysin S and necrotising infections produced by group G streptococcus.
        Lancet. 2002; 359: 124-129
        • Kachlany S.C.
        • Planet P.J.
        • Bhattacharjee M.K.
        • et al.
        Nonspecific adherence by Actinobacillus actinomycetemcomitans requires genes widespread in bacteria and archaea.
        J Bacteriol. 2000; 182: 6169-6176
        • Kurtz S.
        • Phillippy A.
        • Delcher A.L.
        • et al.
        Versatile and open software for comparing large genomes.
        Genome Biol. 2004; 5: R12
        • Leyden J.J.
        The evolving role of Propionibacterium acnes in acne.
        Semin Cutan Med Surg. 2001; 20: 139-143
        • Lomholt H.B.
        • Kilian M.
        Population genetic analysis of Propionibacterium acnes identifies a subpopulation and epidemic clones associated with acne.
        PLoS One. 2010; 5: e12277
        • Makarova K.S.
        • Haft D.H.
        • Barrangou R.
        • et al.
        Evolution and classification of the CRISPR-Cas systems.
        Nat Rev Microbiol. 2011; 9: 467-477
        • McDowell A.
        • Gao A.
        • Barnard E.
        • et al.
        A novel multilocus sequence typing scheme for the opportunistic pathogen Propionibacterium acnes and characterization of type I cell surface-associated antigens.
        Microbiology. 2011; 157: 1990-2003
        • Nizet V.
        • Beall B.
        • Bast D.J.
        • et al.
        Genetic locus for streptolysin S production by group A streptococcus.
        Infect Immun. 2000; 68: 4245-4254
        • Ross J.I.
        • Eady E.A.
        • Cove J.H.
        • et al.
        16S rRNA mutation associated with tetracycline resistance in a gram-positive bacterium.
        Antimicrob Agents Chemother. 1998; 42: 1702-1705
        • Ross J.I.
        • Snelling A.M.
        • Eady E.A.
        • et al.
        Phenotypic and genotypic characterization of antibiotic-resistant Propionibacterium acnes isolated from acne patients attending dermatology clinics in Europe, the U.S.A., Japan and Australia.
        Br J Dermatol. 2001; 144: 339-346
        • Schreiner H.C.
        • Sinatra K.
        • Kaplan J.B.
        • et al.
        Tight-adherence genes of Actinobacillus actinomycetemcomitans are required for virulence in a rat model.
        Proc Natl Acad Sci USA. 2003; 100: 7295-7300
        • Tamura K.
        • Dudley J.
        • Nei M.
        • et al.
        MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0.
        Mol Biol Evol. 2007; 24: 1596-1599
        • Tarr P.I.
        • Gordon C.A.
        • Chandler W.L.
        Shiga-toxin-producing Escherichia coli and haemolytic uraemic syndrome.
        Lancet. 2005; 365: 1073-1086
        • Webster G.F.
        Inflammation in acne vulgaris.
        J Am Acad Dermatol. 1995; 33: 247-253
        • White G.M.
        Recent findings in the epidemiologic evidence, classification, and subtypes of acne vulgaris.
        J Am Acad Dermatol. 1998; 39: S34-S37
        • Zerbino D.R.
        • Birney E.
        Velvet: algorithms for de novo short read assembly using de Bruijn graphs.
        Genome Res. 2008; 18: 821-829