For analysis by HTS, 6 DNA samples extracted from cutaneous swabs obtained from the skin surface of facial areas (forehead and eyebrows) of patients previously studied by PCR for MCPyV sequences were selected (5). These samples included 1 from an index patient with an MCC on his elbow and 5 from the skin of 5 healthy persons.

For investigation of prevalence by specific nested PCR, 120 skin specimens were similarly obtained from 120 volunteers. The median age of these persons was 48 years (range 19–96 years); 30 persons were 57–96 years of age (median age 71 years). This group of 120 volunteers was composed of 40 patients hospitalized or attending outpatient clinics at the dermatology unit at Montpellier University Hospital for various skin disorders (including 8 patients with MCC; median age 75 years, range 57–86 years), 20 immunocompromised patients without skin lesions (10 patients infected with HIV-1 without skin symptoms and 10 renal transplant recipients receiving immunosuppressive regimens [steroids, mycophenolate mofetil, and calcineurin inhibitors]), and 60 healthy controls.

Respiratory samples tested were composed of 46 bronchoalveolar lavage samples obtained from hospitalized patients in intensive care units with acute respiratory failure of unknown origin and 46 nasopharyngeal aspirates from children in the pediatric emergency unit at Montpellier University Hospital with various respiratory tract disorders. An additional 92 fecal samples were obtained from children hospitalized in the pediatric unit for gastroenteritis.

DNA from all samples was extracted as described (5). For HTS, DNA was amplified by using a bacteriophage ϕ29 polymerase-based rolling circle amplification assay and random primers. The protocol of the QIAGEN REPLI-g Midi Kit (QIAGEN, Courtaboeuf, France) was followed as recommended by the manufacturer.

HTS was performed by using the Illumina HiSeq 2000 apparatus (Illumina Inc., San Diego, CA, USA) at GATC Biotech AG (Konstanz, Germany). Five micrograms of high molecular weight amplified DNA was divided into 200–350-nt fragments to which adaptors were ligated. These adaptors included a nucleotide tag that enabled multiplexing several samples per lane or channel. Sequencing was conducted at a mean depth per sample of 8.9 × 10⁶ paired-end reads of 100 nt (range 7.6–10.3 × 10⁶ reads).

Sequences were first sorted by using a subtractive database comparison procedure. Several assembly programs dedicated to short or medium reads were used to generate contigs: Velvet (www.ebi.ac.uk/~zerbino/velvet), SOAPdenovo (http://soap.genomics.org.cn), and CLC Genomics Workbench (www.clcbio.com) (J. Cheval et al., unpub. data). Comparison of single reads and contigs with known genomic and taxonomic data was performed by using dedicated specialized viral, bacterial, and generalist databases created and maintained at the Institut Pasteur (GenBank viral and bacterial databases). Aforementioned databases were screened by using BLASTN and BLASTX (http://blast.ncbi.nlm.nih.gov/Blast.cgi). We used BLAST software (Paracel, Pasadena, CA, USA) capable of executing searches on multiple nonshared memory processors simultaneously.

The entire sequence of the Institut Pasteur polyomavirus (IPPyV) strain genome was analyzed and annotated by using CLC Genomics Workbench (CLC Bio, Aarhus, Denmark). GenBank reference sequences of other members of the family Polyomaviridae used were JCPyV (NC_001699), BKPyV (NC_001538), KIPyV (NC_009238), WuPyV (NC_009539), MCPyV (NC_010277), SV40 (NC_001669), TSPyV (NC_014361), and LPV (M30540). Protein structures were visualized by using Pymol (Delano Scientific LLC, San Francisco, CA, USA).

Phylogenetic reconstructions were based on separate analyses of nucleotide sequences from viral protein 1 (VP1) and large T antigen (LT). A 974-nt region of monkey B-lymphotropic papovavirus (reference sequence M30540.1 from the VP1 coding sequence) was aligned with corresponding regions from the polyomaviruses available in GenBank. For the LT matrix, a 1,453-nt region (same reference sequence as for VP1) was used for analysis. Sequences were aligned by using SeaView version 4 (6) and the Muscle algorithm (7). Only partial but contiguous parts of each alignment were included in final matrixes because alignment in some regions was not possible. The 742–1715 and 2904–4356 nt regions were included in matrixes. The matrixes were based on monkey B-lymphotropic papovavirus M30540 and African green monkey polyomavirus NC_004763 sequence annotations for VP1 and LT genes, respectively.

Phylogenetic analyses were performed by using a probabilistic (Bayesian) approach implemented in BEAST (8). Matrixes were tested against 88 substitution models by using jModelTest software (9). On the basis of results obtained with jModelTest software, the generalized time reversible substitution model (with invariant sites and a gamma site heterogeneity distribution) was used for analysis. The 3-codon partition model of evolution and the Yule speciation process were also specified as priors. A Markov chain Monte Carlo method was used to set 30,000,000 states to obtain an adequate posterior effective sample size >300. Pertinence of nodes was evaluated by using posterior probabilities. Sequences used in phylogenetic reconstructions were obtained from the National Center for Biotechnology Information (NCBI) database (www.ncbi.nlm.nih.gov/nucleotide). Confidence intervals for proportions were calculated according to the efficient-score method (corrected for continuity) (10) (http://dogsbody.psych.mun.ca/VassarStats).

For sequencing the IPPyV by the Sanger method, 9 primer pairs were designed to amplify the entire genome by reference to the contigs assembled from HTS data acquired in the first phase (see Results). Primer sequences and protocols are available upon request. After the genome was sequenced, we developed a specific nested PCR for detection of IPPyV in samples by using primers based on the IPPyV genome sequence and designed by using PrimerPro 3.4 software (www.changbioscience.com): VP1_354F (5′-ACCATATCAGTAGGATAGGTA-3′) and VP1_354R (5′-TGAATTGTATGGCTACAGTGC-3′) for the outer PCR, and VP1_198F (5′-CACTGGGATAGTTCCTGAGG-3′) and VP1_198R (5′-CCTAATGCTACTACCCTCCCT-3′) for the inner PCR. These primers were designed to avoid amplification of other known human polyomaviruses.

The study was reviewed and approved by the Institut Pasteur Comité de Recherche Clinique and the French Commission Nationale Informatique et Libertés (09.465). Consent was provided by participants for obtaining human samples according to French regulations.

Using 8,052,770 Illumina reads obtained from DNA extracted from the skin surface of the MCC index patient, we assembled the complete genome of MCPyV. We found numerous papillomavirus contigs and contigs covering more than half of the genomes of HPyV6 and HPyV7. Additionally, 14 other contigs were assembled that showed a better homology with LPV (NCBI accession no. M30540, version M30540.1, GI:333282) than with any other virus present at that time in the NCBI database, including other human or animal members of the family Polyomaviridae. On the basis of sequence of 8 of the 14 obtained contigs, which were distributed along the LPV genome, we defined a set of 9 primer pairs encompassing the entire target genome. These primers enabled amplification of the entire genome by PCR and analysis of its sequence of 5,028 nt by using the Sanger method, which confirmed the circular nature of the genome.

Whole genome organization of the IPPyV strain, which exhibits general molecular characteristics of polyomaviruses, is shown in Figure 1. It encodes analogs of small T antigen, LT, and structural proteins VP1, VP2, and VP3 and does not appear to encode an agnoprotein. Pairwise amino acid identity was 100% between IPPyV and HPyV9 proteins, 72%–80% between IPPyV and LPV proteins, and much lower for other known family Polyomaviridae members (Table). Because the nucleotide sequence of IPPyV is nearly identical to that of HPyV9, with a difference of only 2 nt in a noncoding region at nt 4449, it appears that IPPyV should be considered a strain of HPyV9. Its sequence have been was deposited in GenBank (accession no. FR823284).

Reconstructions of VP1 and LT phylogenies on the basis of nucleotide sequences clustered the mammalian polyomaviruses and placed the species Avipolyomavirus in basal position when rooting with the oldest known Avipolyomavirus (Budgerigar fledging virus) (Figure 2). Despite the divergences described below, VP1 and LT of HPyV9 were closely related to those of the monkey B-lymphotropic polyomavirus. Moreover, VP1 and LT phylogenies consistently identify several monophyletic groups among mammalian polyomaviruses (Figure 2). Nevertheless, the incongruence of VP1 and LT signals induce notable differences in the topology of these 2 phylogenies. For instance, LT of bovine polyomavirus is closely related to one of the common ancestors of all other mammalian polyomaviruses, whereas its VP1 is closely related to VP1 of the sea lion polyomavirus. Incongruence between VP1 and T phylogenies has been observed for HPyV6 and HPyV7. Nomenclature described in proposals of the International Committee on Taxonomy of Viruses is shown in Figure 2, even though the species Orthopolyomavirus is not monophyletic and therefore should be considered cautiously.

We compared the secondary structure of VP1 from LPV and HPyV9 because the external capsid protein of Polyomaviridae is known to interact with the cell receptor and because antibodies cross-reacting with LPV VP1 have been detected in a large proportion of humans. Overall amino acid identity was 87.1%. The VP1 monomer consists of antiparallel β-strands folded into a jelly roll β-barrel structure. Three outer loops (BC, DE, and HI) are exposed outside the pentamer core and are most likely recognized by antibodies. Using the crystal structure of SV40 VP1 (3BWQ), we mapped the amino acids that differed between the 2 proteins. Polymorphic residues are present in the 3 VP1 loops, and the BC and HI loops appear more conserved than DE loop, which thus shows the major differences (Figure 3).

We first confirmed by specific nested PCR the presence of HPyV9 in the skin swab specimen of the index patient in which the virus had been identified by HTS. We also identified this virus by nested PCR in a second cutaneous sample from the same index patient case obtained 20 months after the first sample. Because HPyV9 was identified in a patient who had MCC, we analyzed the skin surface of 7 other MCC patients. HPyV9 was detected in 1 patient 80 years of age. The overall prevalence in the MCC group was thus 2/8 (25%, 95% confidence interval 4.4%–64%; p = 0.05).

Because we were interested in possible interhuman transmission of HPyV9, a skin swab specimen from the wife of the first index patient was analyzed and found to be positive for this virus. We sampled 111 skin swab specimens from healthy persons or non-MCC patients who had no known contact with MCC patients and screened them with the same nested PCR. Only 1 healthy 30-year-old person had HPyV9, which demonstrated low prevalence in this control group of 1/111 (0.9%, 95% confidence interval 0.05%–5.6%; p = 0.05). Age range of the MCC population was 57–86 years (median 75 years), and because we were unable to detect HPyV9 among the 30 older controls (age range 57–96 years, median 71 years), this finding suggests that the rate of detection in MCC samples was not biased by older ages of these patients. None of the 92 respiratory and 92 stool specimens was positive for HPyV9.