Skip to main content
  • Research article
  • Open access
  • Published:

Genetic associations and shared environmental effects on the skin microbiome of Korean twins



The skin is the outermost layer of the human body and one of the key sites for host-microbe interactions. Both environmental and host genetic factors influence microbial communities in distinct anatomical niches, but little is known about their interplay in shaping the skin microbiome. Here, we investigate the heritable components of the skin microbiome and their association with host genetic factors.


Based on our analysis of the microbiota from 45 individuals including monozygotic and dizygotic twins aged 26–55 years and their mothers, we found that skin microbial diversity was significantly influenced by age and skin pigmentation. Heritability analysis revealed genetic and shared environmental impacts on the skin microbiome. Furthermore, we observed a strong association between the abundance of Corynebacterium jeikeium and single nucleotide polymorphisms (SNPs) in the host FLG gene related to epidermal barrier function.


This study reveals an intimate association of the human skin microbiome and host genes, and increases our understanding of the role of human genetic factors in establishing a microbial ecosystem on the body surface.


The skin provides a protective barrier against invading microbial pathogens, while also serving as a habitat for a plethora of commensal bacteria. The composition of skin microbiota varies among individuals, but a core set of microbes can be found in all subjects for each anatomical location [13]. Recent advances in high-throughput sequencing technologies have enabled rapid analysis of microbial communities indigenous to specific niches in the human body [4]. There is significant inter-individual variation in the skin microbiota, which is attributable to genetic differences and environmental factors. Colonization of the host body by bacteria occurs immediately after birth, with more complex microbial communities developing in distinct anatomical niches afterwards [5]. Succession of the microbiota within individual hosts is also driven by host-intrinsic factors, such as age, gender, genotype and health status, as well as environmental and lifestyle factors, such as climate, light exposure, and detergent use [6, 7]. Previous studies have quantitatively explored host genetic effects on the gut microbiome composition using animal models and humans [811]. These studies revealed a significant influence of host genetic factors on the microbiota. However, the impact of host gene expression on skin microbial diversity in human subjects has not been fully investigated.

The majority of complex human traits are governed by multiple gene interactions [12]. These sophisticated interactions differ among individuals based on environmental factors. Studies on twins are therefore a useful method for estimating the genetic and environmental effects that a given factor has on human health [13]. Twin heritability and linkage studies segregate the genetic and environmental effects based on the assumption that monozygotic twins are genetically identical, and so their phenotypic differences are of environmental origin. Additionally, the shared environmental effects can be determined from monozygotic (MZ) and dizygotic (DZ) twins, as well as parent-offspring pairs [13]. Our knowledge of individual differences in the skin microbiota and genetic markers associated with its composition remains limited. In this study, we analyzed the skin microbiota of 16 MZ and 8 DZ twins between 26 and 55 years of age and their mothers (n = 16) to identify the heritable components of the skin microbiome and their association with host genetic factors. To this end, we focused our analysis on host genes related to key dermatological conditions, including sebum production, skin humidity, pigmentation, epidermal barrier function, and hair follicle development. Using the microbial composition as quantitative traits, we identified one human single nucleotide polymorphisms (SNPs) strongly associated with the abundance of Corynebacterium jeikeium on the skin.

Results and discussion

Study population and composition of the human skin microbiome

The study population included Korean subjects (n = 45) with eight MZ and four DZ twin pairs, five singletons, and 16 of their mothers (Table 1). Skin traits including skin humidity and pigmentation were measured in 32 subjects. Subjects were categorized as having sufficiently moisturized versus dry skin, and mid- versus light-tone skin color based on threshold values in an arbitrary unit for skin humidity (45 AU) and pigmentation (150 AU), respectively. A total of 860,547 sequencing reads from the V2-V3 regions of the human 16S rRNA genes were generated from 45 skin swab samples, for an average of 19,123 reads per sample. We chose the flexor surface of the right upper arm to minimize variations due to differences in the use of cosmetics, surface temperature, acidity, humidity, handedness, and other external factors. The skin microbiota was first assessed at the phylum level; Actinobacteria (50.0 %), Firmicutes (22.0 %), Proteobacteria (16.1 %), and Bacteroidetes (6.0 %) were the dominant phyla among the samples collected. The composition of the skin microbiome in our study was consistent with previous reports [7]. Further classification of the skin microbiota at the genus level revealed high inter-individual variation with Propionibacterium (37.11 %), Staphylococcus (6.79 %), and Streptococcus (6.95 %) prevailing on the skin (Fig. 1a). At the species level, P. acnes was the most abundant inhabitant, accounting for 23.62 % of the skin microbiota. Occupancy by the remaining species accounted for less than 1 % each.

Table. 1 Summary of the study population (N = 45)
Fig. 1
figure 1

Diversity of the skin microbiota and functional traits of the microbial communities. a Taxonomic classification of the skin microbiome at the genus level. Relative bacterial abundance for each individual is shown. A familial relationship is indicated among those individuals represented by the same number. Twin pairs are presented in the parenthesis. b Skin metagenomes predicted using PICRUSt (Methods). Relative contribution of each functional pathway is determined for the collective microbiota of all subjects. Error bars, SEM

To explore the skin microbe-host interaction further, we used PICRUSt [14], an algorithm that detects the functional capabilities of a community by comparing its metagenome with reference genomes using microbial communities identified from all subjects. PICRUSt classified 200 functional pathways from the 16S rRNA sequences of human skin microbiota.

Functional traits determined based on 16S rRNA genes showed categorical similarity with a recent skin metagenomic study (Fig. 1b), with a large proportion of the functions associated with carbohydrate, amino acid, vitamin, and nucleotide metabolism [15]. Discrepancies in the relative contribution of each pathway between the studies could be due to use of different sampling sites and analytical methods. The previous metagenomic study reflects the microbial functions from two male subjects at mid-twenties with samples from five different body parts. Additionally, 16S rRNA based predictions of the functional traits by PICRUSt could possibly cause the discrepancies. Detailed categories of the functional traits and Nearest Sequenced Taxon Index (NTSI) values are provided in Additional file 1: Table S1 and Additional file 1: Table S2, respectively.

Diversity of the human skin microbiome

Microbial diversity “within” and “between” subject groups is denoted by alpha and beta diversity, respectively. In an assessment of species richness, the rarefaction curves, reflective of alpha diversity, showed significantly higher levels of bacterial diversity at ages greater than 35 years (Wilcoxon rank-sum test, W = 109, P < 0.001; Fig. 2). The median value was used as the cutoff age at 35 to evenly disperse individuals into two groups. Previous studies have shown that age is a primary contributor to increased microbial diversity in the skin. A study comparing the skin microbiota between children and adolescents revealed microbial shifts with the onset and progression of sexual maturation [16]. The amount of sebaceous wax ester is known to be closely associated with age, increasing between the ages of 15–35 years and decreasing afterwards [17]. An increase in skin lipids is thought to create an acidic environment, reducing the numbers of acid-susceptible bacteria—such as Staphylococcal and Streptococcal species—on adolescent skin. Disparities in sweat, sebum, and hormone production, as well as lifestyle choices such as the use of cosmetics, may underlie gender differences in skin microbiota [2, 3, 15]. However, gender effects were not present on the inner wrist (Wilcoxon rank-sum test, W = 147, P = 0.144) when mothers were excluded to avoid confounding by age. The gender effects on bacterial diversity seem to depend on the skin site; greater diversity of bacteria in females was evident on hands [2] but not at the base of the neck [18]. Additionally, we could not observe the significant effect of skin humidity (Wilcoxon rank-sum test, W = 126, P = 0.534) on the skin. This is likely due to the relatively low variation in humidity in the skin of the upper arm where the microbiota was sampled. Interestingly, skin pigmentation significantly influenced skin bacteria, increasing the diversity in individuals with an intermediate skin tone (Wilcoxon rank-sum test, W = 176, P = 0.045; Additional file 1: Figure S1). While there have been reports on the effects of bacterial species on melanogenesis [1921], the exact nature of the relationship between the skin microbiota and pigmentation remains to be determined. Our data revealed only a weak association with skin humidity (Pearson correlation, r = 0.393, P = 0.029; Additional file 1: Table S3). To compare the beta diversity profiles across different status of age, gender, skin humidity, and pigmentation, we performed a principal coordinate analysis (PCoA) using weighted and unweighted UniFrac metrics (Additional file 1: Figure S2), wherein the distance between microbial communities was calculated based on phylogenetic information. The results were not indicative of complete separation of microbiotas based on the host phenotypes (ANOSIM, P > 0.05 for all host traits) indicating that the groups divided by age, gender, skin humidity, and pigmentation did not include significantly different microbial compositions.

Fig. 2
figure 2

Differences in skin microbial diversity analyzed using rarefaction curves a Age. b Gender. c Skin humidity (CV, Corneomerty Value). d Pigmentation (MI, Melanin Index). Error bars indicate 95 % confidence intervals. * Significance is based on the Wilcoxon rank-sum test

Heritability of the skin microbiota

Comparison of the skin microbiota between twins and their mothers was performed using weighted and unweighted UniFrac distances (Fig. 3a and b). With both metrics, MZ twins showed the highest similarity, followed by DZ twins and mother-twins. The dissimilarity of unrelated subject was the highest under the unweighted UniFrac metrics whereas it was the second highest next to mother-twins under the weighted UniFrac distance. The similarity of the skin microbiota of MZ twins was comparable to that of DZ twins, but significantly higher than those in unrelated subjects and their mothers (P < 0.05). Our subsequent analysis of the skin microbiota focused on heritability and the additional effects stemming from common environments (Table 2). Heritability ranged from 40.9–56.4 % with the abundance of Roseomonas in the skin showed the highest degree of heritability. Although the role is not well-defined, a recent study comparing skin microbes on the face and other body parts found an increased abundance of Roseomonas on the scapula and inner elbow [22]. Two of the taxonomic branches described in this study were found to possess a heritable component: Corynebacteriaceae/ Corynebacterium and Brevibacteriaceae/ Brevibacterium/ Brevibacterium aureum. Corynebacterium is a major constituent of the normal skin flora along with Staphylococcus, Propionium, Streptococcus, and Pseudomonas [23]. C. jeikeium, the only species belonging to genus Corynebacterium described in this study, is commonly found on skin epithelium, and it has been shown to cause nosocomial infections in immunocompromised patients [24, 25]. Brevibacterium sp. is a Gram-negative bacterium capable of degrading cyfluthrin, an organic compound with xenobiotic characteristics and toxicity to the reproductive, neural, and respiratory systems in humans [26]. The range of shared-environmental effects on the skin microbiome was 56–72.9 %, with Leuconostocaceae/ Weissella/ Weissella cibaria (70.1 %; P = 0.021) and Propionibacteriaceae/ Propionibacterium/ Propionibacterium acnes (57.7 %; P = 0.027) significantly influenced by the common environment from family to species. W. cibaria is a lactic acid bacterium belonging to the family of Leuconostocaceae [27]. Lactic acid bacteria require specific nutrients, such as carbohydrates, amino acids, vitamins, purines, and pyrimidines for growth [28]. Despite these requirements, the bacteria have been isolated from a variety of environments, including the oral cavity, intestinal tract, and fermented foods [28, 29]. This adaptability likely accounts for the presence of the bacteria in two seemingly incongruous locations. P. acnes, the most dominant and well-known skin commensal bacterium [30], is heavily reliant on the metabolism of fatty acids found in sebum secretions. Soap, detergents, astringents, and other factors that alter the amount of sebum may account for the environmental effects on P. acnes colonization [31]. Individual-specific effects ranged from 29.9–59.1 %.

Fig. 3
figure 3

Comparison of the skin microbiota of twins and their families in terms of heritability and environmental effects. a Weighted and b unweighted UniFrac distance between twin pairs (4 UniFrac distance values for DZ pairs and 8 for MZ pairs), twins and mothers (25 UniFrac distance values), and unrelated individuals (909 UniFrac distance values). The error bars represent the standard error. *,P< 0.05 based on Wilcoxon rank sum test

Table. 2 Heritability and household effects of the skin microbiota after adjustment for age and gender

Genetic association of skin microbiotas

We next explored associations between the skin microbiota composition and host genetic variation. For this investigation, we focused on SNPs in a panel of select host genes known to affect sebum production (MC5R, FA2H, DGAT1, DGAT2, SCD1, and ELOVL4), pigmentation (DCT, OCA2, and TYRP1), skin humidity (AQP3), skin barrier function (KLF4, FLG, POF1B, and SPINK5), and hair follicle development (EDA1, EDAR, EDARADD, and PKP1) (Additional file 1: Table S4). The selection of genes was performed manually. One of 275 SNPs was significantly associated with the abundance of a specific bacterial species (Table 3 and Fig. 4a), where significance was defined using a cut-off level of P < 0.05/275 = 0.000182. The identified SNP was localized in a gene known to play a role in skin barrier function (FLG). Marker rs6996321 (allele T) in FLG was negatively associated with C. jeikeium, although the abundance difference was small between the host genotypes (P < 0.05; Fig. 4b). FLG defects are known to cause allergic skin diseases, such as atopic eczema and ichthyosis vulgaris [32], both of which are characterized by dry skin. Correspondingly, possession of the minor allele T tended lower skin humidity but it was not statistically proved (Fig. 4c, P > 0.05). FLG, which encodes filaggrin, is a structural protein in the cornified envelop of the stratum corneum critical for skin barrier function [33]. Mutations in the serine protease St14 in mice, which impairs filaggrin processing, also results in an ichthyotic phenotype [34]. Importantly, Corynebacterium is overrepresented in the skin microbiota of St14-deficient mice, suggestive of a possible link between bacteria and filaggrin processing and pathogenesis. The allelic association with bacterial abundance was further confirmed using LEfSe analysis, an algorithm based on calculation of significantly different phenotypes using the Kruskal-Willis test followed by the Wilcoxon rank-sum test for post hoc analysis (Fig. 4d) [35]. Linear Discriminant Analysis (LDA) scores allowed us to estimate the effect size of the differentially abundant features. LEfSe analysis showed an additional association between Peptoniphilus asaccharolyticus and the minor allele of FLG. Although not identified in the SNP analysis, a substantial percentage of Peptoniphilus abundance could be attributed to heritability, indicating a possible genetic relationship between the bacteria and host. Interactions between the host and the skin microbiota have been analyzed in association with multiple skin diseases [3638]. For example, a recent quantitative trait locus (QTL) mapping of mouse skin microbiotas revealed a strong connection between disease susceptibility and the host genotype-dependent microbial composition [38]. This study, performed on a genome-wide scale, detected QTL regions related to innate immune activation. This and our investigation did not identify overlapping QTL regions, which is conceivable because our SNP analysis was based on the specifically selected gene categories with known implications for human skin structure and function. Beyond the identification of cis-regulatory element associated SNPs, recent SNP association studies are expanding to include flanking regions of the target genes and long-range enhancer/promoter [3941]. In this way, it is able to discover the interplay between target SNPs and their transcriptional regulators and ultimately, practical roles of the SNPs. Further studies on such long-range SNPs should allow us deeper insights into the genetic associations on the skin microbiome.

Table. 3 Associations between the skin microbiota and SNPs of targeted human genes after adjustment for age and gender
Fig. 4
figure 4

Effect of the host genotype on skin traits. a Manhattan plot summarizing the results of an association analysis of 275 candidate SNPs with Corynebacterium jeikeium. Each dot represents the candidate SNP plotted across the genome. b Relative abundance of C. jeikeium with respect to rs6996321. The means and standard errors are indicated by solid and dashed bars, respectively. c Level of skin humidity by genotype at rs6996321. Boxes represent the 25th percentile, median, and 75th percentile. Whiskers represent the lowest values and the highest values of skin humidity. Filled circles represent outliers. AU: arbitrary unit. d LDA (linear discriminant analysis) plot of skin bacteria found by LEfSe showing their association with the host genotype at rs6996321. Cladogram on the right indicates the phylogenetic distribution of the skin bacteria. Each color represents host genotype: CC genotype in red and CT in green. Circles are arranged by phylogenetic levels from kingdom, phylum, class, family, genus, and species from inside out


In this study, we describe host genetic factors and other host-intrinsic characteristics that directly influence the composition of the skin microbiota. We further demonstrated a strong association between the composition of skin microbiota and human genetic factors related to skin barrier function. Analysis of genetically identical MZ twins and half-identical DZ twins can be used to identify and discriminate genetic and environmental effects on the microbiome [42]. Investigation of the twin pairs proved that the skin microbiome is shaped both by host genetics and their environmental factors. Additionally, analysis of human genetic traits associated with the skin microbiota was performed using a candidate-gene approach. For this investigation, we chose five representative traits that affect dermatologic conditions; namely, sebum production, pigmentation, skin humidity, skin barrier function, and hair follicle formation. To confirm the results of the association analysis, the bacterial abundance associated with host alleles was further analyzed using LEfSe. One SNP, located in a gene related to skin barrier function, was significantly associated with C. jeikeium. C. jeikeium abundance was lower in subjects containing the minor allele of FLG, which encodes filaggrin, a structural protein in the cornified envelop of the stratum corneum critical for skin barrier function.

Although we have explored intrinsic characteristics of the human skin microbiome in aspect of host genetics and environmental factors, our study has statistical limitations from a small sample size. In genetic association studies, a sufficient number of samples are critical to detect causality between genes and phenotypes. Furthermore, the collective size of bacteria identified by microbiome analysis creates more stringent p-values to achieve asking for increased sample size to identify additional links between host genetic factors and the composition of the skin microbiota. To overcome such limitations, we tried to assess the genetic impacts put on the skin microbiome using different analytical approaches and the results are supportive of our findings. Yet, more expansive genome-wide association studies with increased sample size are warranted to identify additional links between host genetic factors and the composition of the skin microbiota. Our results provide insight into skin traits associated with individual’s microbiome and potential genomic and microbial targets for skin healthcare.


Study population and sample collection

Study subjects were either MZ or DZ twins, and their mothers were recruited for the Healthy Twin Study as part of the Korean Genome Epidemiology Study [43] between September 2010 and August 2011. Zygosity of twins was confirmed using either 16 short tandem repeat (STR) markers (15 autosomal markers and 1 sex-determining marker) (67 %) or a self-administered zygosity questionnaire (33 %), which showed >90 % accuracy [44]. A total of 51 trios were initially selected, but 6 subjects were excluded due to intake of antibiotics within 3 months of sample collection. Thus, the final sample size was 45 individuals including 8 MZ twin pairs, 4 DZ twin pairs, and 21 family members; the members were composed of 5 parent-offspring pairs and 11 mothers of the twin pairs. Age of the mothers ranged from 49–79 years, and twin children and singletons ranged from 26–55 years. Two subjects exhibited mild symptoms of atopy, but were not excluded as they were only reported to have such symptoms and not prescribed any medication. Dermatological phenotypes such as pigmentation and humidity were available from 32 participants. The participants provided their written informed consent to participate in this study. All experiments involving human subjects were approved by the Korea Centers for Disease Control and the Institutional Review Board of the Seoul National University (IRB No. 144–2011–07–11).

The inner wrist of the right arm was swabbed with two sterile cotton swabs moistened with ST solution (0.15 M NaCl with 0.1 % Tween 20) [45]. The number of strokes was 4–5 times per sampling with gentle pressure. The heads of the cotton swabs were stored in 100 μl of ST solution at −80 °C until use. Swab samples were collected after 3-hr medical checkup securing the time from contact with water and soap. Additionally, skin phenotype—including pigmentation and humidity—were determined using different probes with a C + K multi probe adapter MPA 9 (Courage + Khazaka Electronic GmbH, Köln, Germany) on the same date. Sample collection and visiting was done at set times strictly. Skin pigmentation was measured from the flexor surface of right arm using a Mexameter® MX 18, as described previously [46]. The device quantifies the ratio of light emitted and reflected by skin chromophores and calculates the amount of melanin. Melanin levels were expressed as the melanin index, which ranged between 0 and 150 arbitrary units (AU) for light skin tone, and 150 and 250 AU for mid-tone, skin pigmentation. For skin humidity, capacitance measurement was performed using a Corneometer®CM825 [47]. Each phenotype was measured three times at 18–23 °C and 40–60 % relative humidity. Corneometry values of greater than 45 AU indicate sufficiently moisturized skin, and values less than 45 AU were considered dry skin. The categories of each trait followed the manufacturer’s instructions. The average value of the repeated measurements was used unless the measurements differed by more than 10 %. In this case, the median of the three measurements was taken.

DNA extraction and 454 pyrosequencing

Total genomic DNA was extracted following the bead-beating extraction protocol [48]. Briefly, the cotton swab and ST solution were added to 500 μL of extraction buffer (200 mM NaCl, 200 mM Tris, and 20 mM EDTA; pH 8), 500 μL of phenol:chloroform:isoamyl-alcohol (25:24:1; pH 7.9) (Sigma, Steinheim, Germany), 210 μL of 20 % SDS, and 500 μL of zirconia-silica beads (0.1 mm in diameter; Biospec Products Inc., Bartlesville, OK). The mixture was homogenized using a Vortex Adaptor (MoBio Laboratories, Solana Beach, CA) for 2 min at room temperature. DNA extraction was performed with 500 μL of phenol: chloroform: isoamyl-alcohol (25:24:1; pH 7.9), followed by isopropanol precipitation. The nucleic acid solutions were stored at −70 °C until use. The V2 and V3 regions of the 16S rRNA genes were amplified from total DNA, extracted, and pyrosequenced as described previously [49]. The experiments were quality controlled using proper negative and positive controls at the stage of DNA extraction and PCR amplification, respectively. The cotton swabs without swabbing the inner arm were used as the negative controls. PCR amplification of these controls did not amplify any bacterial DNA, thus they were not used for subsequent sequence analysis. The positive controls used were DNA from E.coli at the stage of PCR amplification. The sequence data have been submitted to the EMBL databases under accession number PRJEB5864 (

Bioinformatic analysis using 16S rRNA sequences

Sequence data were analyzed using the QIIME software package (version 1.5.0) [50]. Before the sequence analysis, quality filtering was performed including removal of low-quality sequences (<200 bp) and ambiguous reads and end-trimming. Subsequently, homopolymers were removed by denoising the sequence data set in the QIIME pipeline [51]. Representative sequence sets were chosen using UCLUST and clustered at the 97 % sequence similarity level. Processed sequences were aligned using PyNAST [52], and taxonomy was assigned using the ribosomal database project (RDP) classifier [53], and the Greengenes Database (gg_97_otus_4feb2011.fasta) was used as the reference [54]. The minimum confidence score for the taxonomy assignment to sequences was 0.8. Chimera sequences (10.55 %) were excluded from downstream analyses prior to the generation of phylogenic trees or OTU tables using ChimeraSlayer algorithm [55]. Bacterial diversity both within and between samples was assessed through alpha using Chao1 measure [56] and beta diversity using weighted and unweighted UniFrac distances [57]. Phylogenetic tree was generated using the FastTree method [58]. UniFrac metrics were also used to identify differences between sample pairs. Analysis of similarity (ANOSIM) was performed to test the statistical significance in the differences. Functional traits were determined from 16S-rRNA-based sequences using PICRUSt-0.9.1( [14]. Unclassified pathways were excluded from the functional analysis.

Heritability analysis

Heritability estimates of each skin microbe were calculated by variance component methods using Sequential Oligogenic Linkage Analysis Routines (SOLAR, version 6.6.2; Southwest Foundation for Biomedical Research, San Antonio, TX, USA) [59]. Skin bacterial abundances at all taxonomic levels were used as quantitative traits. As the bacterial abundances were not normally distributed, inverse normal transformation was applied to the traits before the heritability analysis. Additionally, the monozygotic twins were separately categorized from the rest of participants. SOLAR uses a maximum-likelihood method, which allows incorporation of fixed covariate effects (age, gender, and interaction between age and gender). Also, the software allows working with extended families and complicated pedigree of different age and sex. Variance component analysis decomposes the total variation (σp2) into additive genetic effects (σa2) and residual non-genetic variance (individual effect; σe2), which can be further specified to unmeasured common environmental effects (σc2). We fitted the variance component model into additive genetic and non-genetic components and added the common environmental variance when the effects (σc2) were present (σp2 = σa2 + σc2 + σe2). Thus, the heritability was estimated as the ratio of variance attributed to the additive genetic components and the total variance (σa2 /σp2).

Skin traits associated with SNPs

For SNP genotyping, venous blood from each subject was collected and genomic DNA was extracted using the i-genomic Clinic DNA Extraction Kit (Intron, Seongnam, South Korea). Genotyping was performed on Affymetrix Genome Wide Human SNP Array 6.0 (Affymetrix, Inc., Santa Clara, CA). Exclusion criteria included SNPs with Minor Allele Frequency < 0.01, Hardy Weinberg Equilibrium < 0.001, low call rate (<95 %), Mendelian error, and non-Mendelian error. Mendelian and non-Mendelian errors were identified using the PEDSTATS 0.6.12 [60] and Merlin 1.1.2 [61] software.

Association analysis of skin microbiotas was performed using Plink 1.07 ( [62], as described previously [63]. Briefly, the QFAM (family-based test of association for quantitative traits) model was used to analyze SNPs obtained from skin-trait related genes (Additional file 1: Table S4). Candidate genes were selected based on their roles in regulatory functions or skin metabolism. The list of SNPs was obtained from the Single Nucleotide Polymorphism Database (dbSNP; Build 137) [64]. Skin bacterial abundances at all taxonomic levels were used as quantitative traits. QFAM uses a simple linear regression and corrects for family structure based on adaptive permutation; in this study, 100,000 permutations were performed. This adaptive permutation also alleviates the assumption about normality of data set. Prior to analysis, the bacterial abundances were adjusted for age and gender by fitting to a regression model in R version 3.0.2 [65], as QFAM does not allow for covariates. For MZ twins, only one individual from each pair was analyzed. The alleles of the significant SNPs were examined with the bacterial abundances using LEfSe (linear discriminant analysis [LDA] coupled with effect size measurements) software [35]. The linear discriminant analysis (LDA) effect size greater than 2 was used as the threshold for discriminative bacteria.

Statistical analysis

Species richness and UniFrac distances were analyzed by the Wilcoxon rank sum test (two-tailed) using R version 3.0.2 [65]. Correlation test of skin pigmentation was performed by the Pearson correlation with the SPSS software, ver. 21 (Armonk, NY, US). Association of host genotype with skin humidity was tested using the Kruskal-Willis test.


  1. Staudinger T, Pipal A, Redl B. Molecular analysis of the prevalent microbiota of human male and female forehead skin compared to forearm skin and the influence of make-up. J Appl Microbiol. 2011;110(6):1381–9.

    Article  CAS  PubMed  Google Scholar 

  2. Fierer N, Hamady M, Lauber CL, Knight R. The influence of sex, handedness, and washing on the diversity of hand surface bacteria. Proc Natl Acad Sci U S A. 2008;105(46):17994–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  3. Grice EA, Kong HH, Conlan S, Deming CB, Davis J, Young AC, et al. Topographical and temporal diversity of the human skin microbiome. Science. 2009;324(5931):1190–2.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  4. Kong HH. Skin microbiome: genomics-based insights into the diversity and role of skin microbes. Trends Mol Med. 2011;17(6):320–8.

    Article  PubMed Central  PubMed  Google Scholar 

  5. Dominguez-Bello MG, Costello EK, Contreras M, Magris M, Hidalgo G, Fierer N, et al. Delivery mode shapes the acquisition and structure of the initial microbiota across multiple body habitats in newborns. Proc Natl Acad Sci U S A. 2010;107(26):11971–5.

    Article  PubMed Central  PubMed  Google Scholar 

  6. Roth RR, James WD. Microbial ecology of the skin. Annu Rev Microbiol. 1988;42:441–64.

    Article  CAS  PubMed  Google Scholar 

  7. Gao Z, Tseng CH, Pei ZH, Blaser MJ. Molecular analysis of human forearm superficial skin bacterial biota. Proc Natl Acad Sci U S A. 2007;104(8):2927–32.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  8. Benson AK, Kelly SA, Legge R, Ma F, Low SJ, Kim J, et al. Individuality in gut microbiota composition is a complex polygenic trait shaped by multiple environmental and host genetic factors. Proc Natl Acad Sci U S A. 2010;107(44):18933–8.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  9. Kovacs A, Ben-Jacob N, Tayem H, Halperin E, Iraqi FA, Gophna U. Genotype is a stronger determinant than sex of the mouse gut microbiota. Microb Ecol. 2011;61(2):423–8.

    Article  PubMed  Google Scholar 

  10. Zhao L, Wang G, Siegel P, He C, Wang H, Zhao W, et al. Quantitative genetic background of the host influences gut microbiomes in chickens. Sci Rep. 2013;3:1163.

    PubMed Central  PubMed  Google Scholar 

  11. Goodrich JK, Waters JL, Poole AC, Sutter JL, Koren O, Blekhman R, et al. Human genetics shape the gut microbiome. Cell. 2014;159(4):789–99.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  12. Yang J, Lee T, Kim J, Cho MC, Han BG, Lee JY, et al. Ubiquitous polygenicity of human complex traits: genome-wide analysis of 49 traits in Koreans. PLoS Genet. 2013;9(3):e1003355.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  13. Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley RE, et al. A core gut microbiome in obese and lean twins. Nature. 2009;457(7228):480–4.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  14. Langille MG, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol. 2013;31(9):814–21.

    Article  CAS  PubMed  Google Scholar 

  15. Mathieu A, Delmont TO, Vogel TM, Robe P, Nalin R, Simonet P. Life on human surfaces: skin metagenomics. PLoS One. 2013;8(6):e65288.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  16. Oh J, Conlan S, Polley EC, Segre JA, Kong HH. Shifts in human skin and nares microbiota of healthy children and adults. Genome Med. 2012;4(10):77.

    Article  PubMed Central  PubMed  Google Scholar 

  17. Jacobsen E, Billings JK, Frantz RA, Kinney CK, Stewart ME, Downing DT. Age-related changes in sebaceous wax ester secretion rates in men and women. J Invest Dermatol. 1985;85(5):483–5.

    Article  CAS  PubMed  Google Scholar 

  18. Hyman RW, Babrzadeh F, Palm C, Wang C, Fukushima M, Davis RW. Changes in the human skin microbiome over one year’s time. Am J Microbiol. 2012;3(2):18.

    Article  Google Scholar 

  19. Liu WS, Kuan YD, Chiu KH, Wang WK, Chang FH, Liu CH, et al. The extract of Rhodobacter sphaeroides inhibits melanogenesis through the MEK/ERK signaling pathway. Mar Drugs. 2013;11(6):1899–908.

    Article  PubMed Central  PubMed  Google Scholar 

  20. Huang HC, Chang TM. Antioxidative properties and inhibitory effect of Bifidobacterium adolescentis on melanogenesis. World J Microbiol Biotechnol. 2012;28(9):2903–12.

    Article  PubMed  Google Scholar 

  21. Chen YM, Shih T-W, Chiu CP, Pan T-M, Tsai T-Y. Effects of lactic acid bacteria-fermented soy milk on melanogenesis in B16F0 melanocytes. J Funct Foods. 2013;5(1):395–405.

    Article  CAS  Google Scholar 

  22. Hillion M, Mijouin L, Jaouen T, Barreau M, Meunier P, Lefeuvre L, et al. Comparative study of normal and sensitive skin aerobic bacterial populations. Microbiologyopen. 2013;2(6):953–61.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  23. Cogen A, Nizet V, Gallo R. Skin microbiota: a source of disease or defence? Br J Dermatol. 2008;158(3):442–55.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Coyle MB, Lipsky B. Coryneform bacteria in infectious diseases: clinical and laboratory aspects. Clin Microbiol Rev. 1990;3(3):227–46.

    PubMed Central  CAS  PubMed  Google Scholar 

  25. Larson E, McGinley K, Leyden J, Cooley M, Talbot GH. Skin colonization with antibiotic-resistant (JK group) and antibiotic-sensitive lipophilic diphtheroids in hospitalized and normal adults. J Infect Dis. 1986;153(4):701–6.

    Article  CAS  PubMed  Google Scholar 

  26. Chen S, Dong YH, Chang C, Deng Y, Zhang XF, Zhong G, et al. Characterization of a novel cyfluthrin-degrading bacterial strain Brevibacterium aureum and its biochemical degradation pathway. Bioresour Technol. 2013;132:16–23.

    Article  CAS  PubMed  Google Scholar 

  27. Srionnual S, Yanagida F, Lin L-H, Hsiao K-N, Chen Y-s. Weissellicin 110, a Newly Discovered Bacteriocin from Weissella cibaria 110, Isolated from Plaa-Som, a Fermented Fish Product from Thailand. Appl Environ Microbiol. 2007;73(7):2247–50.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Aguirre M, Collins M. Lactic acid bacteria and human clinical infection. J Appl Bacteriol. 1993;75(2):95–107.

    Article  CAS  PubMed  Google Scholar 

  29. Ahrne S, Nobaek S, Jeppsson B, Adlerberth I, Wold AE, Molin G. The normal Lactobacillus flora of healthy human rectal and oral mucosa. J Appl Microbiol. 1998;85(1):88–94.

    Article  CAS  PubMed  Google Scholar 

  30. Bruggemann H, Henne A, Hoster F, Liesegang H, Wiezer A, Strittmatter A, et al. The complete genome sequence of Propionibacterium acnes, a commensal of human skin. Science. 2004;305(5684):671–3.

    Article  PubMed  Google Scholar 

  31. Leyden JJ. Therapy for acne vulgaris. N Engl J Med. 1997;336(16):1156.

    Article  CAS  PubMed  Google Scholar 

  32. Irvine AD, McLean WI, Leung DY. Filaggrin mutations associated with skin and allergic diseases. N Engl J Med. 2011;365(14):1315–27.

    Article  CAS  PubMed  Google Scholar 

  33. McAleer MA, Irvine AD. The multifunctional role of filaggrin in allergic skin disease. J Allergy Clin Immunol. 2013;131(2):280–91.

    Article  CAS  PubMed  Google Scholar 

  34. Scharschmidt TC, List K, Grice EA, Szabo R, Renaud G, Lee CC, et al. Matriptase-deficient mice exhibit ichthyotic skin with a selective shift in skin microbiota. J Invest Dermatol. 2009;129(10):2435–42.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  35. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12(6):R60.

    Article  PubMed Central  PubMed  Google Scholar 

  36. Fry L, Baker BS, Powles AV, Fahlen A, Engstrand L. Is chronic plaque psoriasis triggered by microbiota in the skin? Br J Dermatol. 2013;169(1):47–52.

    Article  CAS  PubMed  Google Scholar 

  37. Kong HH, Oh J, Deming C, Conlan S, Grice EA, Beatson MA, et al. Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis. Genome Res. 2012;22(5):850–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  38. Srinivas G, Moller S, Wang J, Kunzel S, Zillikens D, Baines JF, et al. Genome-wide mapping of gene-microbiota interactions in susceptibility to autoimmune skin blistering. Nat Commun. 2013;4:2462.

    Article  PubMed Central  PubMed  Google Scholar 

  39. Smemo S, Tena JJ, Kim KH, Gamazon ER, Sakabe NJ, Gomez-Marin C, et al. Obesity-associated variants within FTO form long-range functional connections with IRX3. Nature. 2014;507(7492):371–5.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  40. He H, Li W, Liyanarachchi S, Srinivas M, Wang Y, Akagi K, et al. Multiple functional variants in long-range enhancer elements contribute to the risk of SNP rs965513 in thyroid cancer. Proc Natl Acad Sci U S A. 2015;112(19):6128–33.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  41. Lemire M, Zaidi SH, Ban M, Ge B, Aissi D, Germain M, et al. Long-range epigenetic regulation is conferred by genetic variation located at thousands of independent loci. Nat Commun. 2015;6:6326.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  42. Martinl N, Boomsma D, Machin G. A twin-pronged attack on complex traits. Nat Genet. 1997;17(4):387–92.

    Article  Google Scholar 

  43. Sung J, Cho S-I, Lee K, Ha M, Choi E-Y, Choi J-S, et al. Healthy Twin: a twin-family study of Korea protocols and current status. Twin Res Hum Genet. 2006;9(6):844–8.

    Article  PubMed  Google Scholar 

  44. Song YM, Lee D-H, Lee MK, Lee K, Lee HJ, Hong EJ, et al. Validity of the zygosity questionnaire and characteristics of zygosity-misdiagnosed twin pairs in the Healthy Twin Study of Korea. Twin Res Hum Genet. 2010;13(03):223–30.

    Article  PubMed  Google Scholar 

  45. Paulino LC, Tseng C-H, Strober BE, Blaser MJ. Molecular analysis of fungal microbiota in samples from healthy human skin and psoriatic lesions. J Clin Microbiol. 2006;44(8):2933–41.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  46. Piérard G. EEMCO guidance for the assessment of skin colour. J Eur Acad Dermatol Venereol. 1998;10(1):1–11.

    PubMed  Google Scholar 

  47. Berardesca E. EEMCO guidance for the assessment of stratum corneum hydration: electrical methods. Skin Res Technol. 1997;3(2):126–32.

    Article  Google Scholar 

  48. Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley RE, et al. A core gut microbiome in obese and lean twins. Nature. 2008;457(7228):480–4.

    Article  PubMed Central  PubMed  Google Scholar 

  49. Lee JE, Lee S, Lee H, Song Y-M, Lee K, Han MJ, et al. Association of the vaginal microbiota with human papillomavirus infection in a Korean twin cohort. PLoS One. 2013;8(5):e63514.

    Article  PubMed Central  PubMed  Google Scholar 

  50. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7(5):335–6.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  51. Quince C, Lanzén A, Curtis TP, Davenport RJ, Hall N, Head IM, et al. Accurate determination of microbial diversity from 454 pyrosequencing data. Nat Methods. 2009;6(9):639–41.

    Article  CAS  PubMed  Google Scholar 

  52. DeSantis T, Hugenholtz P, Keller K, Brodie E, Larsen N, Piceno Y, et al. NAST: a multiple sequence alignment server for comparative analysis of 16S rRNA genes. Nucleic Acids Res. 2006;34 suppl 2:W394–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  53. Cole JR, Wang Q, Cardenas E, Fish J, Chai B, Farris RJ, et al. The Ribosomal Database Project: improved alignments and new tools for rRNA analysis. Nucleic Acids Res. 2009;37 suppl 1:D141–5.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  54. Soergel DA, Dey N, Knight R, Brenner SE. Selection of primers for optimal taxonomic classification of environmental 16S rRNA gene sequences. ISME J. 2012;6(7):1440–4.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  55. Haas BJ, Gevers D, Earl AM, Feldgarden M, Ward DV, Giannoukos G, et al. Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons. Genome Res. 2011;21(3):494–504.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  56. Chao A: Nonparametric estimation of the number of classes in a population. Scand J Stat 1984;11(3):265–270

  57. Lozupone C, Knight R. UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol. 2005;71(12):8228–35.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  58. Price MN, Dehal PS, Arkin AP. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS One. 2010;5(3):e9490.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998;62(5):1198–211.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  60. Wigginton JE, Abecasis GR. PEDSTATS: descriptive statistics, graphics and quality assessment for gene mapping data. Bioinformatics. 2005;21(16):3445–7.

    Article  CAS  PubMed  Google Scholar 

  61. Abecasis GR, Cherny SS, Cookson WO, Cardon LR. Merlin—rapid analysis of dense genetic maps using sparse gene flow trees. Nat Genet. 2001;30(1):97–101.

    Article  PubMed  Google Scholar 

  62. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  63. Loukola A, Wedenoja J, Keskitalo-Vuokko K, Broms U, Korhonen T, Ripatti S, et al. Genome-wide association study on detailed profiles of smoking behavior and nicotine dependence in a twin sample. Mol Psychiatry. 2014;19(5):615–24.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  64. Wheeler DL, Barrett T, Benson DA, Bryant SH, Canese K, Chetvernin V, et al. Database resources of the national center for biotechnology information. Nucleic Acids Res. 2007;35 suppl 1:D5–D12.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  65. TR Development Core Team. R: A language and environment forstatistical computing. R Foundation for Statistical Computing. Vienna, Austria;2008.

Download references


This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2010–0029113).

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Joohon Sung or GwangPyo Ko.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

G.K. and J.S. conceived and designed the research. J.Y.S. performed the experiments. J.Y.S. and S.L. analyzed the data, J.Y.S., J.P. and G.K. coordinated the manuscript writing. All authors read and approved the final manuscript.

Additional file

Additional file 1: Figure S1.

Number of OTUs by skin pigmentation. Boxes represent the 25th percentile, median, and 75th percentile. Whiskers represent the lowest values and the highest values of skin pigmentation. Filled circles represent outliers. MI: Melanin Index. Figure S2. Differences in skin microbial diversity. Principal coordinate analysis (PCoA) for (A) age, (B) gender (C) skin humidity (CV, Corneomerty value), and (D) pigmentation (MI, Melanin Index) based on 16S rRNA genes (weighted/unweighted UniFrac distances). Left: weighted, Right: unweighted measures. Table S1. Functional traits and their abundance based on the 16S rRNA genes of members of the skin microbiota. Table S2. Nearest Sequenced Taxon Index (NSTI) values for individuals. Table S3. Correlation between skin pigmentation and selected variables. Table S4. Summary of the SNPs of 14 targeted genes (DOCX 871 kb).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Si, J., Lee, S., Park, J.M. et al. Genetic associations and shared environmental effects on the skin microbiome of Korean twins. BMC Genomics 16, 992 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: