Significance

Skin forms a critical protective barrier between a mammal and its external environment. Baseline data on the mammalian skin microbiome elucidates which microorganisms are found on healthy skin and provides insight into mammalian evolutionary history. To our knowledge, this study represents the largest existing mammalian skin microbiome survey. Our findings demonstrate that human skin is distinct, not only from other Primates, but from all 10 mammalian orders sampled. Identifying significant similarities between branching of mammalian phylogenetic trees and relatedness trees for their corresponding microbial communities raises the possibility that mammals have experienced coevolution between skin microbiota and their corresponding host species.

Abstract

Skin is the largest organ of the body and represents the primary physical barrier between mammals and their external environment, yet the factors that govern skin microbial community composition among mammals are poorly understood. The objective of this research was to generate a skin microbiota baseline for members of the class Mammalia, testing the effects of host species, geographic location, body region, and biological sex. Skin from the back, torso, and inner thighs of 177 nonhuman mammals was sampled, representing individuals from 38 species and 10 mammalian orders. Animals were sampled from farms, zoos, households, and the wild. The DNA extracts from all skin swabs were amplified by PCR and sequenced, targeting the V3-V4 regions of bacterial and archaeal 16S rRNA genes. Previously published skin microbiome data from 20 human participants, sampled and sequenced using an identical protocol to the nonhuman mammals, were included to make this a comprehensive analysis. Human skin microbial communities were distinct and significantly less diverse than all other sampled mammalian orders. The factor most strongly associated with microbial community data for all samples was whether the host was a human. Within nonhuman samples, host taxonomic order was the most significant factor influencing skin microbiota, followed by the geographic location of the habitat. By comparing the congruence between host phylogeny and microbial community dendrograms, we observed that Artiodactyla (even-toed ungulates) and Perissodactyla (odd-toed ungulates) had significant congruence, providing evidence of phylosymbiosis between skin microbial communities and their hosts.

Skin is the largest organ of the body and represents the primary physical barrier between mammals and their external environment. Characterization of skin microbiota is essential for diagnosing skin conditions (1), understanding how an animal coevolves with its microbiota (2), and elucidating interactions between microbiota and the host immune system (3). Skin microorganisms also produce compounds that influence animal physiology, such as intraspecific behavior modifying pheromones (4) and volatile organic compounds that contribute to body odor (5⇓–7).

Although many studies have characterized the human skin microbiome, as reviewed by Grice and Segre (8), far less is known about the skin microbiome of nonhuman mammals, particularly from studies that employed high-throughput sequencing techniques. Multiple studies have been conducted on both wild and captive animals to elucidate the roles that host species, geographic location, body region, and captivity status exhibit on the skin microbiota (9⇓⇓⇓⇓–14). Previous studies of dogs, cats, and bovines have demonstrated that healthy skin microbiota differ between healthy and diseased or allergic animals (13⇓⇓⇓–17). Analyzing skin samples of 63 individuals from five primate species revealed that human axillae were associated with distinct microbiota, with lower overall microbial diversity (9). The authors suggested that differences were due to both human hygiene and host–microbe evolution. Skin biopsies and sloughed skin from free-swimming humpback whales (Megaptera novaeangliae) from the North Pacific, South Pacific, and North Atlantic oceans demonstrated that core genera were present despite large geographical distances (11). However, skin microbial communities also exhibited shifts between geographic locations and whale satiation states throughout their migration. A survey of bats determined that the host species, geographic region, and site were significant factors influencing skin communities (10). Microbial diversity of skin and pouch samples from Tasmanian devils identified a strong influence of geographic location and revealed significant differences between wild and captive populations (18). Together, these studies indicate that both phylogeny and habitat can impact skin microbial communities.

Studies that focus on a wide range of species are important for characterizing microbiota–host eco-evolutionary patterns, or “phylosymbiosis.” Phylosymbiosis is the process by which the phylogeny of host species parallels the ecological relatedness of corresponding microbial communities in a given anatomical location (19). Although phylosymbiosis does not necessarily imply coevolution of hosts and their microbiota, coevolution may be one mechanism underlying observations of phylosymbiosis. Using 16S rRNA gene sequencing for high-throughput microbial community analysis, phylosymbiosis was observed by comparing insect host phylogeny with gut microbial community dendrograms (19). In contrast, demonstration of coevolution between microbial taxa and host species would require additional strain-level resolution offered by alternative gene markers with higher phylogenetic resolution, such as DNA gyrase (20). Using such an approach, coevolution was identified within the mammalian gut (20). No study has yet investigated phylosymbiosis of mammalian skin microbiota, which would be a first step toward demonstrating evolutionary histories between mammals and their skin microbiota.

The objectives of this research were to generate a skin microbiome dataset for the class Mammalia and to identify correlations of skin microbiota with species, geographical location, hygiene, and body region. With respect to species–microbiota relationships, we hypothesized that phylosymbiosis occurs between mammals and their skin microbiota. We tested this by comparing mammalian phylogenies with dendrograms generated from their corresponding skin microbial communities. Humans were removed from the phylosymbiosis analysis because our data reveal stark differences between the skin microbiota of humans and all other mammalian species sampled.

Results and Discussion

Comparison of Human and Nonhuman Skin Microbial Communities.

Bacterial and archaeal microbial communities of 589 mammalian skin samples were characterized by sequencing the V3-V4 regions of PCR-amplified 16S rRNA genes (SI Appendix, Table S1). Altogether, a total of 22,728 unique operational taxonomic units (OTUs) were obtained that corresponded to 44 prokaryotic phyla (Dataset S1). An indicator species analysis determined that all human samples had relatively high proportions of Staphylococcus epidermidis, Corynebacterium, and Propionibacterium acnes (Table 1), which is in agreement with previous studies (8, 21, 22). In contrast, nonhuman animals (“animals”) were associated with soil-associated bacteria, such as Arthrobacter and Sphingomonas, albeit at lower average abundance than the human indicators. This finding was corroborated by a core OTU analysis (SI Appendix, Fig. S1), in which a core OTU was defined as one that was present in a minimum of 90% of nonrarefied samples. All mammalian clades shared six core OTUs, including Arthrobacter, Sphingomonas, and Agrobacterium. Five mammalian orders were analyzed further that contained multiple host species and were not composed of cats, dogs, or humans. Each of the orders, except Perissodactyla, had core OTUs that were not shared with any of the other mammalian orders (SI Appendix, Fig. S1). These core OTUs represent microbiota that persist despite different geographical locations and enclosures.

The presence of a large proportion of soil microorganisms (SI Appendix, Fig. S2) may represent bona fide skin microbial communities of these animals, yet likely also reflects frequent contact of the skin with the external environment. Although the sampling protocol for terrestrial mammals did not include a step to rinse off environmental microbiota, as has been adopted for amphibian microbiota studies (23), future studies might test alternative sampling protocols to access the mammalian skin microbiota to minimize sampling of allochthonous microorganisms. Resident and transient microbiota on skin has been discussed previously for human clinical applications (24, 25). Medically, microorganisms are transient if they can easily be disinfected with antiseptics, or removed with soap and water (25). The sampling methods used in this study would detect both the resident and transient microbiota.

Human samples possessed a distinct microbial community from all other nonhuman mammals, except for several domestic pets from the order Carnivora (Fig. 1). In addition, human skin was significantly less diverse than all other mammalian orders, according to both the number of distinct OTUs (Fig. 2A) and Shannon indices (6.54 vs. 3.96, P < 0.001) (Fig. 2B), which confirms observations from a previous study investigating the microbiota of five primate species (9). Other orders in which microbial communities group tightly together include Diprotodontia (kangaroos), Chiroptera (bats), Rodentia (squirrels), and nonhuman Primates (Fig. 1). A subsequent permutational multivariate ANOVA (PERMANOVA) demonstrated that the factor most strongly associated with community variation for all samples was whether the host was a human (F1,587 = 37.8, P < 0.001) (Fig. 3). Because humans have undergone recent evolutionary divergence from other nonhuman primates, such as orangutans (12–15 Mya) and baboons (21–25 Mya) (26), these results suggest that modern human practices, such as spending the majority of time indoors, frequent bathing, and wearing clothing, may have impacted the diversity and composition of measured skin microbiota. A portion of the higher diversity in nonhuman mammals may be due to an increase in the number of transient microbiota. However, a study of previously uncontacted Amerindians demonstrated that changes in lifestyle, such as living outdoors, resulted in higher diversity (27), lending further support to this finding that modern human practices may be rapidly changing the skin microbiota.

Ordination (PCoA) generated by using the Bray–Curtis dissimilarity metric for each of the three body locations sampled. Samples are colored according to mammalian order. (Inset) Ordination colored according to human and nonhuman samples.

Heatmap summarizing the significant metadata factors correlating with the observed skin microbiota for sampled individuals from mammalian orders. Categories with higher PERMANOVA F statistics have higher variation in community dissimilarity. Gray regions of the heatmap represent categories that do not apply. Samples and categories are clustered according to Bray–Curtis dissimilarities.

By analyzing all samples together, random forest modeling identified that human and animal samples could be distinguished correctly 98.5 ± 1.2% of the time. The OTUs that contributed most to the model include Corynebacterium (2.0%), P. acnes (1.2%), Moraxellaceae (1.2%), and Macrococcus (0.8%). These organisms were all within the top 10 most abundant OTUs in a dataset of all samples. A single female human back sample grouped with the majority of the animal samples because of elevated abundances of Luteimonas, Planomicrobium, and Planococcaceae. The animals that were incorrectly classified were house pets, which had elevated levels of Corynebacterium and P. acnes. The specific house pets that grouped predominately with humans lived exclusively indoors. When all pets were removed from the dataset, humans could be distinguished from animals 99.8% of the time, which is 78.2-fold better than expected by chance.

The effects of mammalian taxonomy, body region, and location were analyzed to elucidate whether these factors contributed to the detected skin microbiota. Mammalian order had the strongest association with the observed variation among animal skin communities (PERMANOVA; F9,502 = 11.3, P < 0.001) (Fig. 3), which was followed by geographic location (PERMANOVA; F4,507 = 9.6, P < 0.001) (Fig. 3). Random forest modeling was also conducted on a dataset of only animals to determine how well intrinsic factors (e.g., host taxonomy) and extrinsic factors (e.g., location) could be classified. Animals could be classified best according to their taxonomic order. This model was correct at classifying animals into their corresponding order 87.8 ± 5.0% of the time, which is 5.9-fold better than expected by chance. Lower taxonomic orders, such as family (86.1 ± 3.9%), genus (84.4 ± 4.7%), and species (83.4 ± 6.7%), were progressively classified less accurately. This weaker classification ability may in part be due to smaller number of samples per group for training the model.

The ability to classify accurately from specific locations may in part be due to the soil that is present in a given habitat. Indeed, a study that analyzed the similarity of skin bacterial communities between salamanders and their environment noted that certain taxa were shared between the skin microbiota and the abiotic environment (28), in part due to contact with forest litter. A previous study noted that the host was the most important factor influencing the skin microbiota of amphibians, whereas geographic location was the second-most important factor (29), which aligns closely with both the PERMANOVA and random forest model findings from this study.

Other factors that have been demonstrated to influence the human skin microbiota, such as individuality, biological sex, and body region, exhibited less of an effect on animals sampled in this study. Both taxonomic order and geographic location were classified more accurately than biological sex (65.2 ± 4.5%), body region (39.9 ± 6.0%), or individual animal (36.6 ± 36.7%). This inability to classify individuals is in contrast to human studies that have shown that individuality is one of the most important factors influencing the human skin microbiota, although several of these studies used >15 samples per individual (30⇓–32). It is therefore still possible that the skin microbiota are relatively unique among individual animals, but this cannot be observed with only three samples per animal.

To address whether biological sex influenced the skin microbiota within a species, cat (n = 48), dog (n = 35), and horse (n = 68) samples were analyzed because they contained a relatively large number of sampled individuals and a similar number of samples from both biological sexes. Biological sex was not a significant factor for any of these species (PERMANOVA; cat: F1,15 = 1.15, P = 0.20; dog: F1,11 = 0.79, P = 0.77; horse: F1,20 = 0.94, P = 0.44). All of the domestic cats and dogs were neutered or spayed, so this would not have an influence on the analysis. Although the horses had both neutered/spayed and intact individuals, whether the horse had been neutered or not did not exhibit a significant difference on the skin microbiota (PERMANOVA; F1,20 = 1.09, P = 0.34). The only animal where biological sex explained the most variation was the red kangaroo (PERMANOVA; F1,16 = 2.21, P = 0.002), which was also analyzed because this species exhibited a visual separation between males and females in an ordination (principle coordinates analysis, PCoA) (SI Appendix, Fig. S3). Larger variations in the microbiota among different body regions were observed within Perissodactyla (PERMANOVA; F2,121 = 4.26, P < 0.001) (Fig. 3) and Proboscidea (PERMANOVA; F2,12 = 2.38, P = 0.02) (Fig. 3) than within other orders. The overall low effect from body region is likely due to the body regions sampled. The back, inner thigh, and torso are all covered with hair. A previous study on dogs demonstrated that fur-covered regions had higher species richness and diversity compared with mucosal surfaces (13). Therefore, sampling mucosal surfaces would be expected to result in more distinct differences between body regions within a species.

To ensure that the importance of the host’s taxonomic order on skin microbiota was not overly influenced by orders with fewer samples and locations, the three orders with a large number of samples—Artiodactyla, Carnivora, and Perissodactyla—were analyzed with all other orders removed. Each of these orders had samples sourced from six to eight different locations. Removing orders with fewer samples increased the influence of order (PERMANOVA; F2,385 = 15.1, P < 0.001) and decreased the impact of geographic location (PERMANOVA; F3,384 = 8.7, P < 0.001). Therefore, the effect on microbial communities exhibited by the mammalian host exists despite varying geographic locations, and cannot be fully attributed to certain species only being sampled in a single location.

Pets and Humans.

Although the majority of animals possessed skin microbial communities that were distinct from humans, a subset of pets grouped with humans in ordination space (Fig. 1). In particular, of the 17 animal samples that grouped with humans, 15 were from indoor housecats, whereas the remaining two samples belonged to the backs of dogs that were frequently bathed and groomed (SI Appendix, Table S2). In total, 75% of these samples belonged to animals that were owned by humans who were also sampled for this study. All cats with similar microbial communities to humans had at least two of the three sampled body locations possess this “human” community composition (SI Appendix, Table S2), whereas the two dogs only possessed the human microbial community on their backs. The remaining 11 dogs had similar communities to the other animals, as did all cats that lived outdoors on farms and a single cat that lived in the city, without a dog (Fig. 1, Inset). Interestingly, 11 of 12 indoor cats that lived with a dog possessed similar microbial communities to the other animals. It may be that owning a dog results in an influx of soil microbiota into the home, which in turn is transferred to indoor cats through contact with the built environment or personal contact. Indeed, previous studies show that owning a dog shifts the human microbiota as well as microorganisms detected on built environment surfaces (33, 34). Because our study only sampled skin from 13 dogs and 19 cats, future studies might include a larger sample size of animals that would help further elucidate how owning a dog impacts the skin microbiome of other inhabitants within a household.

Predicted functions based on detected prokaryotic taxa were determined using Functional Annotation of Prokaryotic Taxa v.1.0 (FAPROTAX) (35), revealing several conserved functions on mammalian skin (Fig. 4). Many of these match functions were predicted from human samples from the Human Microbiome Project consortium that were characterized by metagenomic sequencing (36, 37). Animal symbionts and human pathogens were expected because the samples were derived from mammalian hosts. Urea is a component of sweat and provides a nitrogen source, which could explain ureolysis as a predicted skin function, although not all animals sweat in hair-covered regions (38).

There were several functions that were significantly different between humans and animals. In accordance with lower diversity (Fig. 2), humans had significantly fewer predicted functions (34.2 ± 8.4 compared with 51.8 ± 9.4 in animals; P < 0.001). Animals may also have a higher number of predicted functions because they had more soil bacteria, which may have more annotated predicted functions than skin bacteria. Humans had elevated levels of predicted manganese oxidation. Human sweat contains 100 ppb manganese on average, which would result in ∼200–300 mg of manganese secreted each day (39). This concentration is low compared with other trace metal elements, such as zinc and copper (40). Manganese oxidation was predicted to occur from the core human OTU P. acnes (41). In contrast, animals had higher levels of predicted functions involved in the nitrogen cycle and single-carbon compound degradation. Methanol oxidation was attributed to the core OTUs affiliated with Arthrobacter, and methylotrophy with Methylobacterium OTUs, according to the FAPROTAX database. Nitrogen respiration was associated with numerous organisms, such as Paracoccus and Pseudomonas. Understanding if skin microbiota contribute to nitrogen or manganese cycling is important because deficiencies in manganese and nitrogenous compounds, such as nitric oxide, have been associated with delays in wound healing (42). Metabolites have also been linked to aging in mammals due to changes in oxidative stress (43).

Predicting functions based on taxonomy is the first step to elucidating how biochemical processes from skin microbiota are influencing host skin health. Many advances have been made in gut microbiome research that have demonstrated microorganisms within the gastrointestinal tract can influence digestion (44), provide vitamins and amino acids that the host is unable to synthesize alone (45), and impact neurological function through the gut–brain axis (46). Although it is currently unknown how the biogeochemical processes of skin microbiota influence their host, determining which processes are occurring on skin is a first step toward understanding how these microbial functions impact dermatological health. Future studies using metagenomic sequencing may help confirm which of these predicted conserved microbial functions are core to mammalian skin, versus functions that are variable among different host species.

Phylosymbiosis Is Evident Within the Orders Perissodactyla and Artiodactyla.

Phylosymbiosis would result in closely related clades of animals being associated with similar microbial communities (19). This effect can be measured using Robinson–Foulds and matching-cluster values, which compare the congruence between two relatedness trees. A score of 0 indicates that the trees are identical, whereas a score of 1 indicates there is no congruence between the two trees. The matching-cluster method is statistically more robust than Robinson–Foulds and provides a wider distribution of scores between 0 and 1 because it calculates the total composition of the subtrees (47). Previous studies demonstrate that shifts in microbial communities have matched host evolution within insects (19, 48), which were more apparent at the 99% OTU clustering threshold. In this study, comparisons were made at both the 97% and 99% thresholds, using Bray–Curtis, unweighted, and weighted UniFrac distance metrics.

Comparing host mammalian phylogeny to dendrograms of corresponding skin microbial communities for each mammalian host species indicates that relationships among skin microbiota from animals of the orders Perissodactyla and Artiodactyla were similar to their corresponding host evolutionary relationships (Fig. 5). Perissodactyla exhibited phylosymbiosis with all thresholds and distance measures; the only discrepancy in each test case was the skin microbial communities of domesticated horses and Przewalski’s horses (Fig. 5A). The split between the equestrian and rhinoceros clades cannot be attributed to differences in location, such as farm or zoo habitats, because the Przewalski’s horses were sourced from the Toronto Zoo. Artiodactyla (Fig. 5B) possessed a normalized matching cluster score of 0.38 (Table 2), and demonstrated significant congruence with the host phylogeny and microbial dendrogram with both the Bray–Curtis and weighted UniFrac metrics. The largest discrepancy was noted with the sequences from the goats, which grouped with the giraffe and reindeer rather than the sheep. In addition, the host species did not cluster according to the geographic locations from where they were sourced. In contrast, the order Carnivora (Fig. 5C) did not exhibit significant phylosymbiosis and the cat and dog clades did not have distinct microbial communities for either metric of relatedness or distance metric. This observation did not change when all cat and dog samples were removed (Table 2). Therefore, the microbial dendrograms were not being unduly influenced by household animals that undergo frequent grooming and spend the majority of time indoors. It is possible that phylosymbiosis may be more strongly observed within clades of animals that share similar diets or have similar management. All of the sampled animals within Perissodactyla and Artiodactyla were herbivores that graze on local grasses or hay. In contrast, the animals within the Carnivora order had a more diverse diet among member species, such as herbivorous giant pandas, carnivorous lions, and pets that were fed an omnivorous diet. Similar to how diet influences the gut microbiota (49), it may be that the skin microbiota is impacted by diet, as has been shown for skin microbial communities of amphibians (50).

No significant phylosymbiosis was observed when all mammalian orders and humans were analyzed together, except for the unweighted UniFrac measure analyzed using normalized Robinson–Foulds values at the 99% threshold (Table 2). Although, in this single case, animals could be matched significantly better than 100,000 randomized trees of the 38 species, there was very little congruence observed (Table 2), as indicated by the normalized Robinson–Foulds score of 0.93 and the normalized matching cluster score of 0.99. When humans were removed from this dataset, congruence increased modestly, although the host phylogeny and bacterial dendrogram still exhibited little congruence.

Although previous studies have been able to demonstrate phylosymbiosis, they did so under highly controlled laboratory conditions and with fecal samples (19, 49, 51). Skin represents a more transient environment that is influenced by shedding and contact with other surfaces. The animals in this study had several confounding factors, such as different locations and age (Dataset S2). In our study, it is possible that more distinctive congruence would be observed if mammals were sampled at similar time points in their life history and inhabited the same geographic location. Additionally, the potentially transient soil microorganisms that were more abundant on nonhuman mammalian skin may mask phylosymbiosis when all sequenced OTUs are being considered as a community in the phylosymbiosis analysis (SI Appendix, Fig. S1). Future studies should potentially sample before and after washing the skin to observe how this treatment would influence the analysis. We postulate that reducing the number of transient auxiliary organisms from the environment would strengthen the finding of phylosymbiosis because allochthonous microorganisms that do not coevolve with a host would be removed from the analysis.

Archaea Are Present on Mammalian Skin.

Including all human data, archaea comprised only 6,509 of all 6,550,625 nonrarefied sequences analyzed in this study (0.1%) (Dataset S1). Several archaeal clades were represented, including salt-tolerant Halobacteria, the methanogen Methanobrevibacter, and the ammonia-oxidizers of the Thaumarchaeota (Dataset S1). Methanogens likely represent fecal contamination because Methanobrevibacter spp. are the dominant archaea present in the gut (52). However, Halobacteria and thaumarchaeotes, such as Nitrososphaera, have the potential to be resident skin microbiota. Halobacteria are able to tolerate the salt concentrations from sweat (53). Putative ammonia-oxidizing archaea have also been reported on human skin (54, 55). Even though few archaeal sequences were present in our data overall, archaeal sequence proportions, normalized by the number of individuals sampled, were unevenly distributed among sampled animal species (SI Appendix, Fig. S4). For example, averaged archaeal sequence proportions were highest for cape elands (26.1% of all summed archaeal species averages), olive baboons (12.9%), the sable antelope (10.9%), and bovine (5.8%). The methanogens from the phylum Euryarchaeota were the dominant archaeal clade, which was expected because the animals with the most archaeal sequences were predominately ruminants. Averaged thaumarchaeotal sequences were highest for groundhogs (0.9% of all summed archaeal species averages), the swamp wallaby (0.6%), olive baboons (0.5%), and the pony (0.5%). Together, these results expand the known range of putative ammonia-oxidizing archaea on skin beyond human beings.

Our study likely underestimated archaeal abundance and diversity because of primer mismatches to archaeal 16S rRNA genes (56). Indeed, the performance of the Pro341F/Pro805R primer set was analyzed using TestPrime v1.0 on the SILVA database (SI Appendix, Table S3). Only 64.8% of archaeal 16S rRNA genes had zero mismatches to the primer set used in this study. Only 11.9% of ammonia-oxidizing thaumarcheota reference sequences had zero primer mismatches, in contrast to 85.7% of all bacterial sequences. When the number of mismatches was increased to two, 94.9% of all archaea and 95.5% of thaumarcheotes were targeted. A recent study on the gut microbiota of great apes that used both universal prokaryotic and archaea-specific primers determined that the distribution, diversity, and prevalence of archaea in mammalian gut samples is underestimated by up to 90% (56).

Limitations.

The majority of the animals were collected based on opportunistic availability. For example, animals from the Toronto Zoo were sampled during routine veterinary check-ups. The nature of sample collections resulted in an inability to collect from an equal number of representatives from each host taxonomic order and species. For example, the following species only had a single representative animal sampled: alpaca, beaver, pony, sheep, sable antelope, spotted hyena, swamp wallaby, and two-toed sloth. Although we recognize that no significant conclusions can be made about a single host animal within a species, these animals were included in the analysis to maximize coverage of each mammalian order. Much work remains to be conducted within each species to determine intraspecific effects of individuality, body region, and biological sex.

Animal skin microbiota were sampled throughout the calendar year and all samples were frozen until DNA extraction. It is possible that the skin microbiota of outdoor animals may undergo seasonal shifts, especially between the relatively cold winter and warm humid summer in Canada; however, this cannot be tested using a single sampling time for each animal. Future investigations should sample the same individuals temporally to determine if changes in temperature and resulting skin secretion levels influenced microbial community composition and diversity. Moreover, the significant difference in geographic location that was observed may be more pronounced if animals with greater geographic distance were sampled. All of the animals were sampled in Southern Ontario. Sampling the same species from multiple continents is hypothesized to result in more pronounced variations in communities according to location due to significant changes in extrinsic factors, such as soil microorganisms.

A limitation of any amplicon study is biases that arise from primer selection. This study targeted the V3-V4 16S rRNA region. The relative abundances of common skin organisms, such as Propionibacterium and Staphylococcus, may differ in studies that select another portion of the 16S rRNA gene, such as the V1-V3 region (57). Although these biases have been discussed in detail elsewhere (57⇓⇓–60), the V3-V4 region can accurately represent skin microbial communities (61) and is therefore expected to have influenced the results only minimally.

The rodents collected in this study were sourced from both wild and deceased animals. Although samples from deceased animals still grouped with the remaining live animals (Fig. 1), high skin microbial diversity (Fig. 2) may be in part related to initial changes in skin community from decomposition (62). Nonetheless, deceased rodents were collected within 24 h of death and did not have any visible injuries that would result in internal microorganisms from the gastrointestinal tract contaminating the skin.

Conclusion

Our findings demonstrate that human skin is distinct, not only from other primates, but from all 10 mammalian orders sampled. Human samples were dominated by S. epidermidis, Corynebacterium, and P. acnes. Given the recent evolutionary divergence of humans as distinct species from other nonhuman primates, these results suggest that modern human practices, such as living within a built environment, wearing clothing, and washing with soap, have strongly impacted the diversity and composition of the skin microbiota that can be sampled with sterile swabs. Microbial communities from nonhuman mammals were significantly more diverse and were associated with higher levels of traditional soil-associated OTUs that likely represent a combination of resident skin microbiota and transient organisms collected temporarily from their enclosure or natural habitat. Despite these environmental influences, phylosymbiosis was detected in two herbivorous orders: Perissodactyla and Artiodactyla. These results provide first evidence that skin microbial communities assemble in a species-specific manner on mammalian hosts, which may be initial evidence that coevolution is occurring between mammalian hosts and skin microbiota, presumably through maternal inheritance.

Materials and Methods

Ethics.

This study was approved by the Office of Research Ethics (ORE) at the University of Waterloo (A-15-06). Detailed ethical consent for the human portion of this study was reported previously under ORE no. 20993 (63). The following minimally invasive procedures were conducted in accordance with the approved documentation and no animals were harmed throughout this study.

Sample Collection.

Species from 10 orders of the class Mammalia were sampled to characterize the distribution of microorganisms on skin (SI Appendix, Table S1). Both males and females were included for each species, when available, to account for variations in hormone levels and secretions that are known to affect microbial communities (21). A spectrum of habitats and hygiene practices were also included, ranging from frequent grooming of pets and farm animals, to animals in captivity and the wild. Complete information on the biological sex, age, diet, location, health history, grooming, and antibiotic administration were collected (Dataset S2). The inclusion of animals that had received antibiotics in the previous 6 mo was minimal (32 animals; 18%), of which 10 received antibiotics in the 2 mo preceding sampling (5.6%). The results indicate that this exposure did not influence measured microbial diversity of those individuals within each species. Additionally, animals that had received antibiotics did not group separately within ordination space from those that had no history of antibiotics. Therefore, these animals were included to strengthen analyses, such as phylosymbiosis and PERMANOVA. Animals that received antibiotics include the alpaca, aoudad sheep, bovine, olive baboons, cats, dogs, horses, and the Indian rhinoceros.

Animals were sampled from multiple locations in Southern Ontario from November 2015 to September 2016: The African Lion Safari, Kitchener–Waterloo Humane Society (KWHS), Toronto Zoo, pet owners, and from local farms sourced from the University of Guelph. Animals from the Toronto Zoo and African Lion Safari were sampled when they were brought in for regular husbandry practices. Additional companion animals were obtained from volunteers who were recruited by word of mouth. The KWHS supplied wild animals that were collected by KWHS staff within 24 h of death; the specimens were stored in plastic bags in a −20 °C freezer until sampled.

The back, torso, and inner thigh regions of 177 nonhuman mammals were collected using sterile foam swabs (Puritan) according to a previously published protocol (63). In addition, we included data from 77 equivalent samples (the right and left inner thigh were sampled from each participant) from 20 human participants from a previous study that were sampled from November 2015–February 2016 (63) in the analysis for comparison purposes [Sequence Read Archive BioProject ID PRJNA345497 (https://www.ncbi.nlm.nih.gov/sra/PRJNA345497)], for a total of 589 samples. These regions were chosen to capture both moist and dry regions and avoid sensitive areas that may cause distress. Skin was swabbed by moving aside hair or fur with gloved hands to expose the skin. While applying moderate pressure, the skin was swabbed 10 times in a forward and backward motion. The swab was rotated and repeated in adjacent areas for a total of 40 strokes per swab. When the area was complete, the samples were returned to their initial plastic storage container and frozen at −20 °C until further use. All volunteers and veterinary technicians were trained with a detailed protocol to ensure sample collection consistency.

Sample Preparation.

All DNA extractions, PCR protocols, and Illumina sequencing were conducted according to a previously published protocol (63) to enable comparisons between the human and nonhuman samples. In brief, DNA was extracted using the PowerSoil-htp 96-well DNA Isolation Kit (MO BIO Laboratories) and stored at −20 °C until further use. The V3-V4 region of the 16S rRNA gene was amplified using the Pro-341Fi (5′-CCTACGGGNBGCASCAG-3′) and Pro-805Ri (5′-GACTACNVGGGTATCTAATCC-3′) primers (64). Each amplification was performed in triplicate to minimize potential PCR bias from low biomass samples (65), and was conducted in a PCR hood that was UV-treated for 30 min after having undergone a treatment with UltraClean Lab Cleaner (MO BIO Laboratories) to remove DNA, RNA, DNase, and RNAses (66). No template controls, kit extraction controls, and sterile swab controls were included, in addition to including the same samples across all three MiSeq runs to ensure lane-to-lane consistency (“run controls”). Additional descriptions of control methods and results are included in SI Appendix (SI Appendix, Supplementary Information Text and Fig. S5) and Dataset S3. The Pro341F/Pro805R primer set was analyzed using TestPrime v1.0 against the SILVA database to assess potential sequence mismatches against bacterial and archaeal reference sequences (SI Appendix, Table S3) (67, 68).

Processing of Sequence Data.

Raw DNA sequence reads were processed using the same open-source bioinformatics pipeline described previously (63) that was managed by Automation, eXtension, and Integration Of Microbial Ecology (AXIOME) v1.5 (69). PANDAseq v2.8 (70) generated paired-end sequences using the default parameters of a 0.9 quality threshold, a minimum sequence overlap of 10 bases, and a minimum read length of 300 bases. Quantitative Insights Into Microbial Ecology (QIIME) v1.9.0 (71) was used to analyze sequence data, which underwent de novo clustering and chimera/singleton removal at both 97% and 99% cluster identity using UPARSE (72). PyNAST v1.2.2 (73) was used to align 16S rRNA gene sequences. Subsequently, RDP v8.1 (74) assigned prokaryotic taxonomy using Greengenes database v13.8 (75). Samples were rarefied to 1,654 sequences in the dataset that contained all samples (Dataset S1). Analyses such as α-diversity and PERMANOVA underwent 1,000 iterations of rarefication to avoid underrepresenting diversity. Rarefication plots demonstrated that conducting multiple rarefications to determine diversity prevented a loss in diversity levels. Other mammalian skin microbiome studies have used a similar level of sequences to analyze their data (9, 13, 76, 77).

Statistical Analyses.

The majority of statistical analyses were conducted using the same programs and software version numbers as described previously (63). In brief, α-diversity of all samples was measured with the following QIIME commands: multiple_rarefaction.py, alpha_diversity.py, collate_alpha.py, and compare_alpha_diversity.py. Subsequently, the 42 metadata categories were compared using the Bonferroni correction to avoid false-positives due to testing a high number of hypotheses.

β-Diversity was visualized using ordinations generated with the Bray–Curtis distance metric. Ordinations were created in RStudio (78) with the phyloseq (v1.14.0) and ggplot2 (v2.1.0) packages. β-Diversity was measured using PERMANOVA with the Adonis function from the vegan package (v2.4–0) in R. Using 1,000 permutations, the percent variation explained by each metadata category was calculated and visualized in a heatmap using ggplot2, vegan, Heatplus 2.16.0 (79), and RcolorBrewer v1.1.2 (80).

The functions from the prokaryotic clades were predicted using FAPROTAX (35). This conservative algorithm currently matches 80 functions, such as fermentation and methanogenesis, against 7,600 functional annotations of 4,600 prokaryotic taxa.

Phylosymbiosis Analysis.

The phylosymbiosis analysis of skin microbiota profiles and host phylogeny was adapted from a previously described protocol (19). Downloaded COX1 sequences were aligned with Muscle v3.8.31 (81) and edited by removing gap positions and 5′/3′ end overhangs with Jalview v2.9 (82). The final edited alignment was created using RaxML online Blackbox server v8.2 (83). All mammalian host trees were verified to be in concordance with well-established mammalian phylogenies (84⇓⇓⇓–88).

Microbiota dendrograms were constructed using the QIIME v1.9.0 jackknifed_beta_diversity.py command. Species were rarefied to the highest possible sequence count that included all species within the specific taxonomic ranking. This resulted in a rarefication of 1,900 sequences for all mammals, 9,100 for Artiodactyla, 25,700 for Carnivora, and 37,500 for Perissodactyla. Rarefication was conducted 1,000 times and a consensus tree built to correct for the low number of rarefied sequences. Each of the above mammalian clades had skin microbiota consensus dendrograms created at 97% and 99% OTU identity threshold using Bray–Curtis, weighted and unweighted UniFrac distance metrics.

Congruencies between host phylogenies and skin microbiota dendrograms were quantified by calculating normalized Robinson–Foulds (89) and normalized matching-cluster scores. These scores were calculated by using the ape R package (90) and a custom python script created by Brooks et al. (19), respectively. The significance of these scores was determined by constructing 100,000 randomized trees with identical leaf nodes to the skin microbiota dendrograms and comparing each to the host phylogeny to calculate the number of stochastic dendrograms with equivalent or better Robinson–Foulds or matching cluster scores.

Acknowledgments

We thank Rahgavi Poopalarajah, Mayar Zawawi, and Elena Dybner for assistance with sample kit preparation and data entry; Katja Engel for support with Illumina library preparation; Jonathan Ross for assistance with Python; Charles Gray and his team from the African Lion Safari; Dr. Diego Gomez from the University of Guelph; Andrea Dada, Dawn Mihailovic, and the veterinarian technicians from the Toronto Zoo; Amanda Hawkins and Joe Growden from the Kitchener–Waterloo Humane Society; Heather Cray, Emilie Spasov, and all farmers and pet owners who participated in the study; and Andrew Brooks for comments that improved the phylosymbiosis analysis. A.A.R. was supported by a Canada Graduate Scholarship from the Canadian Institutes for Health Research. This research was supported by a Discovery Grant to J.D.N. from the Natural Sciences and Engineering Research Council of Canada.