Recent advances in molecular technologies have identified a much greater diversity of cutaneous flora than what was revealed using culture-based methods 
. Bacterial microbiome studies were mostly based on 16S ribosomal RNA metagenomics since the 16S rRNA
gene contains both consensus regions targeted by PCR and variable regions allowing direct taxonomic assignment whenever possible and phylogenetic studies if required. In the present study, in order to characterize the skin viral microbiota, HTS, a highly comprehensive method based on random sequencing of the entire DNA present in a given sample was applied to superficial skin samples collected by cutaneous swabbing. This method was previously successfully used for studying environmental media, intestinal microbiota and, more recently, human blood samples 
. However, HTS requires large amounts of DNA and comprehensive reference genome sequences for microbial populations, a prerequisite that might raise some concern regarding its application to skin samples. In our study, difficulties regarding the template DNA amount in skin swab samples were circumvented by a first round of whole genome amplification using multiple displacement amplification (MDA), a technique known for its ability to introduce low bias for linear DNA amplification 
. Nervertheless, this method may favor the amplification of circular molecules which might have led to an overestimated proportion of circular DNA genomes (Polyomaviridae, Papillomaviridae, Circoviridae
) as illustrated by the high proportion of eucaryotic viral sequences detected in the two samples 100066 and 100069. Eventually, sequences analysis was restrained to taxonomic assignment with no attempt to obtain any data regarding functional analysis. Also, the contribution of each species was not corrected by genome sizes: such standardization would have lowered by several orders of magnitude the contribution of bacterial species compared to viral species. Thus, comparison of frequencies should be made within bacteria or viruses.
Taxonomic assignment of bacterial sequences detected in our study confirms previous data derived from 16S rRNA analysis. Indeed, most of the bacterial sequences belonged to four dominant phyla, Actinobacteria, Firmicutes, Bacteroidetes and Proteobacteria. Previous molecular analysis of the skin bacterial microbiota have emphasized that bacterial diversity depends on the body site with additional interpersonal and temporal variability 
. More specifically, surface areas could be divided in dry, moist or sebaceous environments regarding skin physiology and it was demonstrated that these conditions are likely to influence the composition of the bacterial microbiome 
. In our study, inter-individual variations are obvious when considering the relative proportion of bacterial sequences within the 6 skin samples tested (Fig S1
). The samples were all collected from the face of the patients, which can be considered as a sebaceous environment, and Propionibacterium
)-related sequences, theoretically expected in such an environment, were indeed detected. Additionally, our study identified numerous sequences from Staphylococcus
spp and Streptococcus
spp (Firmicutes) and Proteobacteria
phylum that are more likely to be found in dry and moist skin environments respectively 
. This discrepancy may be explained by technical limitations of 16S rRNA studies with variable performances according to the phyla under scope whereas HTS metagenomic randomly targets all genomic material with a higher sensitivity and is probably less liable to such variations. In addition, the domination of the 16S based criteria in bacterial taxonomy compared to other molecular characters to date less represented in databases may also contribute to this observed discrepancy.
Detection of microorganisms on human skin is not restricted to bacteria since viruses, fungi and some acarids like Demodex
spp are present on skin of healthy individuals. Most described fungal organisms identified on healthy skin are related to Malassezia
spp whereas Candida
spp is infrequently detected on normal human skin, contrasting with mucosal membrane where this yeast is physiologically present in limited amounts 
. These concepts are supported by our study in which many reads corresponded to Malasseziales
species whereas Saccharomycetales
related reads were rare.
The presence of archaea on the skin has not been previously reported but they might have been missed by 16S rRNA sequencing. Conversely, our results suggest that archaea sequences can be detected on skin surface although the low number of reads raise concerns about the relevancy of these data that need to be confirmed by further investigations. Alternatively, since superficial skin is permanently exposed to surrounding environments as illustrated by the detection of numerous reads related to sequences previously reported from environmental metagenomic studies, we could not rule out contamination by environmental sequences or DNA present in the reagents. Interestingly, sequences from Escherichia coli and Pseudomonas spp species that are in our experience common contaminants of our pipeline did not contribute significantly to the results, most likely because of the abundance of the skin microbiome.
Among viral sequences, only few reads were related to phages. The relative proportion of the various phage species in the 6 samples is in keeping with the spectrum of their bacterial hosts usually present on healthy skin surface. This is in accordance with the detection of bacteriophages in a recently reported virome study conducted in HIV patients peripheral blood, a result likely to be associated to the presence of circulating bacteria in this setting 
. However, the diversity of the phage sequences identified in our study is of great interest because these viruses are involved in lateral gene transfer that can influence the bacterial diversity but also some of the properties of bacteria including virulence and resistance to antibiotics.
According to our data, eucaryotic DNA viruses detected on skin samples displayed a high diversity including various representatives of Papillomaviridae, Polyomaviridae and Circoviridae. However, as mentioned above, some technical limitations must be pointed out that might have reduced a full appraisal of this diversity. First, the initial DNA enrichment step through phi29 MDA might have favored amplification of circular DNA genome as suggested by the recovery of some plasmid sequence. Furthermore, the procedure of DNA extraction and amplification of skin swabs did not allow the recovery of RNA virus genomes. Eventually, as swab samples were all collected from the same facial area, skin site variations, which are known to influence bacterial microbiota, could not be investigated.
Viral colonization of healthy skin was primarily and extensively described for papillomaviruses and the commensal nature of these viruses is now a widely admitted data. Indeed, the presence of resident HPVs on normal skin was reported for both Beta- and Gammapapillomaviruses 
. Our results confirm these data and provide evidence for an asymptomatic carriage of numerous HPV strains (up to 17) on each individual skin sample. This skin carriage may reflect a viral shedding from cutaneous micro-reservoirs that might be hair bulbs, since such high prevalence and multiplicity was mainly described in forehead hair bulb studies 
. Furthermore, we identified 13 new Gamma-HPV strains including 9 strains for which the whole genome sequence was further confirmed by the Sanger method allowing the description of 6 new HPV species (). This suggests that the diversity of the resident cutaneous gamma-HPV group might be actually larger than previously described 
. Chronic asymptomatic shedding of HPVs from skin surface is a well-known feature and most of these viruses are considered as innocent bystanders with no significant tissular damages resulting from their replication which is mainly supported by keratinocytes differentiation.
A similar chronic carriage and shedding at the surface of healthy skin appears to be a hallmark of human polyomaviruses with cutaneous tropism as well. Such a viral ecosystem was suggested by previous reports on MCPyV, HPyV6 and HPyV7 
but no systematic metagenomic study was carried out with high performance methods on normal-appearing skin. In the present study, a correlation was observed between the amount of reads and the MCPyV viral load (evaluated by qPCR in original samples), as illustrated by the relative coverage rates of assembled whole genomes from the six skin samples under scope (). This observation supports the concept of metagenomic approaches being broader and at least as sensitive as specific target amplification and allowing an exhaustive coverage if the corresponding sequences are present at sufficient levels 
. Of note, these results were obtained with a mean depth per sample of 8.6×106
reads, which is far lower than the depth currently permitted by the improvements of Illumina technology that can reach 108
reads. Whole genome comparison of the six MCPyV strains described in our study revealed high sequence similarities to each others and the strain derived from the healthy skin of the MCC patient did not present any specific characteristic. Furthermore no genomic mutation nor deletion were identified suggesting the likely episomal nature of the viral sequences detected in all cases.
Detection of both HPyV6 and HPyV7 sequences in respectively 3 and 1 sample(s) is in keeping with previous studies that have already detected these viruses on superficial layers of healthy skin 
. Accordingly, both of them appear to be less frequently detected than MCPyV and, unlike MCPyV, no skin disease has been associated with these two new viruses to date 
. As for MCPyV, full genome sequences of HPyV6 and HPyV7 identified in this study only demonstrated minor single nucleotide polymorphisms when compared to published data. This low diversity of HPyVs species that appears more limited in its extent than that of HPVs, needs to be confirmed by the analysis of larger cutaneous sample set.
The presence of numerous Circoviridae
representatives on the superficial layers of the skin can be expected since these viruses are highly ubiquitous. These data are strikingly reminiscent of the situation of the very close Anelloviridae
family including TT virus and TTV-like mini virus that are frequently detected in various body sites 
. In our study, the detected sequences are related to the Cyclovirus
genus. Cycloviruses infect mainly animals with likely cross species transmissions 
. Animal cycloviruses could be also detected in feces of primates including humans 
which suggests that the specificity of the detection of cycloviruses on the skin remains to be established.
The biological role of such a highly polymorphic viral microbiota on skin surface is largely hypothetical and has been much less extensively studied than for its bacterial counterpart. Skin bacterial populations are either transient or symbiotic and are likely to provide protection against more virulent organisms, thus resulting in a function of surrogate immune sentinel. A similar function of skin viral microbiota may be hypothesized but warrants to be further investigated. Indeed, resident viruses may act protectively against other bacterial, fungal or viral pathogens through either direct (phages) or indirect (by maintaining a high level of efficiency of skin immune cells) pathways. Similarly, the influence of disturbances in skin bacterial and fungal populations has been emphasized in skin disorders involving local skin immune system dysregulation including seborrheic dermatitis, acne, atopic dermatitis, or psoriasis 
, and the influence of such variations of the viral flora on skin homeostasis and alterations may also be worth considering. Although most cutaneous HPVs and the more recently identified cutaneous HPyVs can be considered as resident symbiotic organisms on normal skin, some of them are associated with skin diseases like HPV5 or HPV8 or HPV50 involved in epidermodysplasia verruciformis or squamous cell carcinoma 
. Also, numerous HPVs and MCPyV are associated to pre-malignant and malignant skin conditions in immunocompromised as well as in immunocompetent patients 
. The hypothesis of the viral skin flora being affected by the immune status of the host is consistent with our experimental data featuring a higher viral diversity of both HPV and HPyV species in the patient with MCC, a rare malignancy mostly affecting elderly and immunocompromised populations 
. Indeed, despite the limited number of cases under scope precluding any statistical analysis, 17 different HPV strains as well as most of the HPyVs previously suspected to have a cutaneous tropism (MCPyV, HPyV6, HPyV7 and HPyV9 ) were identified only on the MCC patient's skin, suggesting that the relationship between skin conditions and/or immunologic status on one hand and the skin virome composition on the other hand should be further investigated.
In summary, this pioneer study is the very first one that applies highly powerful HTS to the description of the viral skin microbiota. Despite a limited number of patients under scope, it demonstrated an unexpected high diversity of DNA viruses on the normal appearing skin with a large interpersonal variability. These viruses are represented essentially by various species of betapapillomaviruses and gammapapillomaviruses, polyomaviruses and circoviruses. Although viruses are generally considered as pathogen agents, our findings highlight the complex viral flora at the surface of the non pathological skin of all individuals. The dynamics and anatomical variations of the skin viroma and their potential variations according to pathological conditions remain to be further studied. It should be also investigated if these viruses, alone or in combination, could represent potential triggers of skin cell proliferation and oncogenesis.