Genetic and Structural Analysis of Merkel Cell Polyomavirus Large T Antigen from Diverse Biological Samples

In 2008, Feng et al. [1], using digital transcriptome subtraction, discovered a new polyomavirus DNA clonally integrated in about 80% of Merkel cell carcinoma (MCC) cases, and therefore named the virus Merkel cell polyomavirus (MCPyV). MCC is a rare and aggressive neuroectodermal malignancy arising in sun-exposed skin of elderly and immunosuppressed patients [2,3]. Nevertheless, MCPyV DNA was also found in several non-MCC sites, from MCC patients and healthy individuals, especially in skin and oral samples [4,5,6,7], although the role of MCPyV in non-MCC sites is unknown.

Large T antigen (LT-ag) polyomaviruses are multifunctional proteins responsible for gene expression regulation and viral DNA replication, with potential oncogenic activity, as seen for the prototype LT-ag from SV40 and MCPyV [8,9,10]. Sequencing of the integrated MCPyV genome from several MCC and cell lines revealed various point mutations in the LT-ag N-terminal DnaJ domain that preserves the LXCXE motif, and a stop codon before the original binding domain (OBD). These truncated forms of MCPyV LT-ag are unable to bind to the viral ORI and therefore to initiate bidirectional replication, suggesting these mutations disrupt LT-ag ability to replicate the integrated viral genome [10,11,12].

Although the consequences of LT-ag truncation for MCPyV oncogenicity have been broadly documented, the effects of other nonsynonymous mutations for MCPyV biology remain obscure and scarcely explored. Therefore, the objective of this study is to analyze the genetic profile of MCPyV LT-ag complete second exon sequences, and to investigate the possible influence of nonsynonymous mutations in protein function through in silico analysis.

The whole second exon of LT-ag (861-3,080 bp) of eleven oral MCPyV-positive samples (8 of saliva and 3 of oral tissue biopsies) from renal transplant recipients and nontransplanted individuals were sequenced using the BigDye™ terminator chemistry kit (Applied Biosystems, Foster City, Calif., USA) according to previously established conditions [5]. They were analyzed along with the whole second exon of all nontruncated sequences available at the NCBI databank to date, with the aid of Molecular Evolutionary and Genetics Analysis (MEGA 5.2) software. Cloned sequences were not used due to potential bias in the number and distribution of mutations.

Overall, thirty-three LT-ag complete second exons were retrieved, comprised of 16 from skin (7 from healthy skin of non-MCC patients; 3 from healthy skin of MCC patients; 3 from non-MCC skin lesions and 3 from MCC lesions), 3 from urine samples, 4 from lung cancer tissue, 3 from nasal swabs, 3 from intestines and 4 from peripheral blood mononuclear cells (PBMC; table 1). The sequence MCC350 (JN038583) was used as a standard truncated LT-ag, and the HF consensus genome (JF813003) as a nontruncated prototype sequence. This study was approved by the Institutional Committee on Research Ethics (protocol 1138/2005).

The MCPyV LT-ag sequences were extracted from the UniProtKB database (UniProtKB accession No. E2IPT4). However, the numbering scheme within this article also considers the first exon (78 amino acids). Homologs of the protein with a solved 3D structure were researched using the Basic Local Alignment Search Tool (BLAST) and the HHPred server (http://toolkit.tuebingen.mpg.de/hhpred). For the helicase domain, a homology model was built with the program Modeller (http://salilab.org/modeller/) using the structure of the SV40 helicase domain as a template (PDB 1SVL). For the retinoblastoma protein binding site (pRb) domain, a homologue with a solved 3D structure was not found, and therefore a model was built using the QUARK structure prediction automated server (http://zhanglab.ccmb.med.umich.edu/QUARK). For the molecular dynamics study, the OBD Merkel structure with PDB ID 3QFQ was used as a wild-type structure. For the molecular dynamics study, all structures of mutant proteins were built with the Modeller program.

Energy minimization and molecular dynamics simulations were performed with the GROMACS package version 4.5.4 [13], using the amber99sb force field [14]. MD simulations were carried out on 20 cores of 2.4 GHz 32 GB of our Linux cluster. The TIP3P water model [15] was used to solvate the protein in an octahedral water box. All simulations were performed in periodic boundary conditions. Bond stretches involving hydrogen atoms were constrained using an LINCS algorithm [16]. The electrostatic interactions were treated by a combination of the particle mesh Ewald method and a switch function. Simulations of solvent molecules and counterions were performed at 310 K and 1 atm for 100 ps, with the protein atoms restrained with harmonic potentials. The equations of motion were integrated using the Verlet Leapfrog algorithm [17]. Systems were minimized for 2,000 steps of steepest descent and then heated from 50 to 310 K using 6 blocks of calculation in which the temperature was gradually increased. Each block lasted 20 ps, totaling 120 ps of the heating process for each system. The equilibrated systems were simulated using a time step of 0.002 ps. All systems were simulated for 20 ns.

In regard to the oral MCPyV samples obtained in this study, the profile of synonymous mutations found in LT-ag were similar between transplanted versus nontransplanted individuals and saliva versus oral biopsy (fig. 1). Only two nonsynonymous mutations, S133T (1,024 nt) and R204T (1,237 nt), were found in saliva samples from transplanted individuals (TxS 1 and TxS 2). In contrast, the mutation pattern between MCPyV LT-ag Brazilian sequences and those retrieved from the NCBI databank from diverse biological samples revealed few similarities (fig. 1). Also, Brazilian and Asian sequences apparently formed homogeneous polymorphic profiles each, which did not occur in other sequences.

With the exception of the pRb region, nonsynonymous mutations were found in all LT-ag domains and inter-domain areas (table 1), and its frequency was significantly higher in samples from lesional tissues (MCC, psoriasis, Kaposi sarcoma, lung cancer, squamous cell carcinoma) than in healthy tissue (p = 0.0018; data not shown). LT-ag sequences of healthy skin from MCC patients are of particular interest given the high number of nonsynonymous mutations found, equivalent to that of the lesional samples.

In order to gain some insight into the possible roles of these mutations in the protein function, a molecular dynamics analysis of the predicted protein structure was performed. The pRb domain is predicted to be intrinsically disordered. Not surprisingly, a close homologue with a solved 3D structure was not found, either by a sequence search or by protein threading methods, and therefore an ab initio model was built using the QUARK structure prediction automated server. According to the model, all of the mutated residues are located in loops at the protein surface, distant from the LXCXE motif, and also far-off the only predicted alpha helix (Aa 80-90). Four out of ten mutations in this region involve either a proline or a glycine residue (table 1), and would therefore modify the flexibility of the pRb domain. The same high percentage of mutations involving glycine or proline residues was found on the loop regions between domains (table 1).

The structure of the OBD domain is already known (PDB 3QFQ). Five mutations were found in this domain (table 1). Residue A337 is located inside the protein in a hydrophobic core. The A337T mutant stability was studied by molecular dynamics simulations. Its structure and dynamics were not different from that of the wild-type domain. Concerning K351R, due to the similarities between lysine and arginine residues, it is plausible that there is no effect of this mutation on protein function. However, this family of proteins can undergo acetylation at lysine residues, a mechanism important for T-ag stability and DNA replication [18,19]. The mutation N386D found in lung cancer tissue could be expected to impact on protein function. Located at the monomer-monomer interface (by comparison with SV40), the mutation of this residue could affect hexamerization [20].

At the helicase domain, 13 nonsynonymous mutations were found. Due to either the residue location in the protein structure or the similar properties of wild-type and mutant residues, some mutations can be hypothesized to not influence protein function: T527S, M531T, V681A, M728T, T753I and I798V. Four mutations, namely F448S, S597R, F714S and L736F, as with the wild-type domain, were submitted to 20 ns of molecular dynamics simulations and no significant differences in protein stability or structure were observed. A swinging of the N-terminal region was observed in all protein models, produced by a conformational change at residues 525-527. Finally, some samples of healthy skin in patients with MCC lesions presented a cysteine-arginine substitution in residue 749. This mutation is located at the p53 binding interface (by comparison to SV40) [21]. As can be seen in figure 2, an arginine residue at position 749 could make a hydrogen bond with Arg438 in p53, increasing the interaction energy between these proteins.

In this study, the analysis of the second exon of the MCPyV LT-ag from non-MCC oral samples revealed no stop codon or any high impact mutations. This result was expected as truncated LT-ag has been described as a tumor-specific signature found in MCC [1,5,10,12] and in other malignancies [22]. Nevertheless, the position of the synonymous mutations detected in all oral samples suggests a pattern and possibly a phylogenetic proximity, regardless the immunosuppressive status and the sample category. Similar results have been found by other authors when studying different biological samples, where genetic polymorphisms were described in MCPyV DNA from urine, nasal swab, PBMC [5] and malignant blood cells [22] from different subjects or the same subject [23], suggesting a strain-specific signature.

Looking at figure 1, it appears that the MCPyV LT-ag polymorphic profile has a geographical influence, as observed by others [24,25]. However, the role of intra-host selective forces in shaping MCPyV DNA remains to be determined, as well as the possibility that certain strains might influence the course of the infection [25]. It should be considered that our apparent mutation profile may also have been influenced by methodological differences between studies, especially the sampling process, number of samples and the characteristics of the participants. However, our scarce knowledge regarding MCPyV phylogeny and pathogenesis limits further conclusions. Yet, it seems that the MCPyV LT-ag polymorphic profile has a geographical conservation.

We could hypothesize that the higher number of nonsynonymous mutations in MCPyV LT-ag from lesional tissues is the consequence of viral adaptation and perhaps escape from host immunity. However, since all Asian (Japanese) samples were collected from lesional tissues, it is difficult to attribute the frequency of mutations to any pathological event.

The in silico study suggests that the mutations found in the second exon of LT-ag would not significantly affect protein functions, regardless the clinical status of the host and the sample type. Also, no apparent pattern of mutations was related to any specific lesion. Stable DNA viruses, like polyomaviruses, are not subjected to strong adaptive selection in the way that high impact mutations are frequently eliminated from the population in favor to neutral/low-impact mutations [24,26]. Therefore, it is tempting to speculate that MCPyV-associated pathologies are ‘all-or-nothing' processes, relying on rare LT-ag truncation events to take place, although the role of other viral genes such as VP1, ST-ag and the NCCR region in MCPyV pathogenesis deserves investigation.

To the best of our knowledge, this work presents the most detailed investigation of the structural and functional properties of LT-ag nontruncated mutations in MCPyV sequences from diverse biological samples. As the information on MCPyV LT-ag structure (mainly the pRb domain), interaction partners and MCPyV pathogenesis increases, the data presented here will gain new light, and the function of nonsynonymous mutations will be better understood.

The authors acknowledge the financial support of FAPERJ.