Identification of virus-encoded microRNAs in divergent Papillomaviruses.

MicroRNAs (miRNAs) are small RNAs that regulate diverse biological processes including multiple aspects of the host-pathogen interface. Consequently, miRNAs are commonly encoded by viruses that undergo long-term persistent infection. Papillomaviruses (PVs) are capable of undergoing persistent infection, but as yet, no widely-accepted PV-encoded miRNAs have been described. The incomplete understanding of PV-encoded miRNAs is due in part to lack of tractable laboratory models for most PV types. To overcome this, we have developed miRNA Discovery by forced Genome Expression (miDGE), a new wet bench approach to miRNA identification that screens numerous pathogen genomes in parallel. Using miDGE, we screened over 73 different PV genomes for the ability to code for miRNAs. Our results show that most PVs are unlikely to code for miRNAs and we conclusively demonstrate a lack of PV miRNA expression in cancers associated with infections of several high risk HPVs. However, we identified five different high-confidence or highly probable miRNAs encoded by four different PVs (Human PVs 17, 37, 41 and a Fringilla coelebs PV (FcPV1)). Extensive in vitro assays confirm the validity of these miRNAs in cell culture and two FcPV1 miRNAs are further confirmed to be expressed in vivo in a natural host. We show that miRNAs from two PVs (HPV41 & FcPV1) are able to regulate viral transcripts corresponding to the early region of the PV genome. Combined, these findings identify the first canonical PV miRNAs and support that miRNAs of either host or viral origin are important regulators of the PV life cycle

Validation of miDGE via confirmation of miRNAs encoded by Japanese Monkey Herpesvirus (JMHV).

<p><b>(A)</b><i>Top</i>: Coverage plot of DNA-seq of expression fragment library containing a cosmid with ~40kb of the Japanese Monkey Herpesvirus (JMRV) genome. <i>Center and bottom</i>: small RNA coverage plots mapping of small RNA-seq reads from HEK293T cells transfected with the expression library (<i>center</i>) or from JMRV-infected rhesus primary fibroblasts (<i>bottom</i>). <b>(B)</b> Small RNA-seq coverage plot in linear (top panels) or logarithmic (bottom panels) across the miRNA-coding region. The genomic location of the 15 JMRV pre-miRNAs is indicated by black arrows at the top. <b>(C)</b> Detailed depiction of small RNA-seq coverage across the first three pre-miRNAs of the JMRV miRNA cluster.</p

High-scoring PV pre-miRNA predictions.

<p>High-scoring PV pre-miRNA predictions.</p

MiRNA candidates identified from miDGE analysis are detected by Northern blot.

<p><b>(A)</b> Northern blot analysis of total RNA from HEK293T cells transfected with indicated miRNA or putative miRNA expression vectors, with ethidium bromide stained low molecular weight RNA shown as a load control. This figure represents a single membrane that was probed first for control SV40 miRNA, then stripped and re-probed for each of the indicated miRNAs. Solid and outline arrows on left side correspond to pre-miRNAs and mature miRNAs, respectively, and approximate sizes are noted on right side. <b>(B)</b> Northern blot analysis of total RNA from HEK293T cells transfected with anti-Drosha siRNA or negative control (NC) siRNA, then re-transfected after 48 hours with respective siRNAs and indicated putative miRNA expression vectors. Total RNA was harvested after 48 hours. Ethidium bromide stained low molecular weight RNA shown as a load control. This figure represents a single membrane that was probed first for control SV40 miRNA, then stripped and re-probed for each of the indicated miRNAs. Membrane was additionally probed for HSUR (Herpesvirus saimiri U RNA 4) as a transfection control. Normalized to HSUR RNA, all indicated miRNAs are detected at higher intensity in negative-control-treated cells than Drosha knockdown cells, suggesting canonical dependence on Drosha processing. The numbers below each blot correspond to the ratio of mature miRNA signal in control cells compared to mature miRNA signal in Drosha-knockdown cells after each signal has been normalized to load control HSUR signal. Note: due to sequence similarity between putative HPV-derived miRNAs, cross-reactivity is seen in hpv17-miR-H1 lanes and hpv37-miR-H1 lanes. <b>(C)</b> Northern blot analysis of RNA from either HEK293T or Dicer KO (NoDice [<a href="http://www.plospathogens.org/article/info:doi/10.1371/journal.ppat.1007156#ppat.1007156.ref068" target="_blank">68</a>]) cells transfected with the indicated miRNA expression vectors. Solid and outline arrows on left side correspond to pre-miRNAs and mature miRNAs, respectively. The numbers below each blot correspond to the ratio of mature miRNA signal to pre-miRNA signal, normalized to the same ratio in HEK293T cells. This figure represents a single membrane that was probed, then stripped and reprobed with the indicated oligonucleotide probes (SV40, FcPV1 F1 and F2, HPV41 miRNA probes) and a second membrane for HPV17 and HPV37 miRNAs. Ethidium bromide stained low molecular weight RNA is shown as a load control.</p

FcPV1 miRNAs are detected <i>in vivo</i> during viral infection.

<p><b>(A)</b> DNA-seq coverage for FcPV1 genome in miDGE library is plotted on log10 scale on the y-axis. The x-axis corresponds to position in the genomic sequence. A schematic of the FcPV1 gene organization taken from NCBI reference sequence NC_004068 is provided at the top. <b>(B)</b> Small RNA-seq coverage for FcPV1 genome in miDGE library is plotted on log10 scale on the y-axis. Values above the x-axis correspond to the forward strand and those below correspond to the reverse strand. The x-axis corresponds to position in the genomic sequence. Peaks corresponding to the newly identified miRNA genes are labeled. <b>(C)</b> Plot of read coverage and start sites for reads mapping to the FcPV1 genome in library prepared with RNA from infected chaffinch leg lesion samples. On the y-axis, coverage is plotted with gray lines and read start counts are plotted with black impulses. Values above the x-axis represent the forward strand and those below represent reads mapping to the negative strand. Genomic position is indicated on the X-axis. Peaks corresponding to the newly identified miRNA genes (fcpv1-miRs-F1 & F2) are labeled as well as the lower expressed candidate miRNA “fcpv1-miR-F3”. The inset photograph is of one of the chaffinches with characteristic leg lesions used in preparation of the small RNA-seq libraries. <b>(D)</b> Plot of read coverage and start sites for reads mapping to the FcPV1 genome in library prepared with RNA from chaffinch pectoral muscle samples. On the y-axis, coverage is plotted with gray lines and read start counts are plotted with black impulses. Values above the x-axis represent the forward strand and those below represent reads mapping to the negative strand. Genomic position is indicated on the X-axis. Peaks corresponding to the newly identified miRNA genes are labeled.</p

miDGE identifies PV-encoded miRNA candidates.

<p><b>(A)</b> Neighbor-joining tree calculated on alignment of L1 nucleotide sequences of the better-covered papillomaviruses (n = 73 with >95% DNA coverage) included in miDGE library. Color indicates genus membership, with miRNA-encoding papillomaviruses in bold and select high risk cancer-associated papillomaviruses in italics. <b>(B)</b> Small RNA coverage distribution of top-scoring miDGE miRNA candidates that were predicted by MiRDeep2. <b>(C)</b> Structures of PV pre-miRNAs were predicted by minimal free energy folding using the RNAfold algorithm. The positions of mature miRNAs observed in small RNA-seq libraries are indicated in red. <b>(D)</b> The position of identified viral pre-miRNAs is denoted by the hairpin shape. The identified seed sequence matches are noted at their respective positions with the sequences of the miRNA and potential targets.</p

PV-encoded miRNAs are active in RISC.

<p>RISC reporter assays for the PV-encoded miRNAs where HEK293 cells were co-transfected with a firefly luciferase transfection control and <i>Renilla</i> luciferase reporter with either perfectly complementary sequence matches for each indicated miRNA, or its respective negative control seed complement mutant. Either of these were co-transfected with control empty miRNA expression vector (blue), the relevant PV miRNA-expression vector (orange), or negative control irrelevant miRNA expression vector (SV40) (gray). Average <i>Renilla</i> luciferase activity relative to firefly luciferase normalized to empty miRNA expression vector control is shown for fcpv1-miR-F1 (N = 3), fcpv1-miR-F2 (N = 5), hpv41-miR-H1 (N = 5), hpv17-miR-H1 (N = 3), and hpv37-miR-H1 (N = 4). Statistical test performed was a Two-Sample t Test. The average <i>Renilla</i> luciferase activity normalized to firefly luciferase activity is shown, error bars indicate Standard Error, and asterisks indicate statistical significance, (*) p≤0.05; (**) p≤0.01.</p

PV-encoded miRNAs can regulate transcript sequences in early genes.

<p>RISC reporter assays for the PV-encoded miRNAs. In all panels, the average <i>Renilla</i> luciferase activity normalized to firefly luciferase activity is shown, error bars indicate Standard Error, and asterisks indicate statistical significance, (*) p≤0.05; (**) p≤0.01. <b>(A)</b> HEK293T cells were co-transfected with either a control empty miRNA expression vector (control) or the indicated PV miRNA expression vector along with both the normalization control firefly luciferase vector and a <i>Renilla</i> luciferase-based reporter plasmids with vector UTR (Empty 3’UTR) or FcPV1 genomic DNA containing both putative miRNA docking sites (FcPV1 Early), N = 4. Statistical test performed was a One-Sample t Test. <b>(B)</b> HEK293T cells were co-transfected with either the SV40 miRNA expression vector (Control) or the indicated PV miRNA expression vector and both the control firefly luciferase vector and the <i>Renilla</i> luciferase-based reporter plasmids with vector UTR (Empty), FcPV1 E1 genomic sequence (E1), or the seed sequence mutant (E1 mut), N = 4. Statistical test performed was a One-Sample t Test. <b>(C)</b> HEK293T cells were co-transfected with either the SV40 miRNA expression vector (Control) or the indicated PV miRNA expression vector and both the control firefly luciferase vector and the <i>Renilla</i> luciferase-based reporter plasmids with vector UTR (Empty), FcPV1 E2 genomic sequence (E2), or the seed sequence mutant (E2 mut), N = 12. Statistical test performed was a One-Sample t Test. <b>(D)</b> HEK293T cells were co-transfected with either the SV40 miRNA expression vector (Control) or the HPV41 miRNA expression vector and both the control firefly luciferase vector and the <i>Renilla</i> luciferase-based reporter plasmids with vector UTR (Empty), HPV41 genomic DNA containing both putative miRNA sites (HPV41 Early), the site in E1 (HPV41 E1), the seed sequence mutant (E1 mut) N = 7 or the site in E2 (HPV41 E2), N = 4. Statistical test performed was a One-Sample t Test.</p

JMRV pre-miRNA predictions.

<p>JMRV pre-miRNA predictions.</p

High Risk HPVs do not encode canonical microRNAs.

<p><b>(A)</b> Coverage plots of miDGE DNA-seq and small RNA-seq across the genomes of HPV18 and HPV17. Filled green plots at the top of each panel show DNA-seq coverage, the three plots underneath show mapped small RNA-seq from: <i>PV</i>: HEK293T-cells transfected with our papillomavirus miDGE library, <i>PV filtered</i>: same reads as in PV, but filtered to eliminate low-complexity reads <i>JMRV</i>: Serving as a negative control, derived from 293T-cells transfected with our JMRV miDGE library. JMRV read counts were normalized to correct for different sequencing depths between PV and JMRV miDGE experiments (see total read counts in <a href="http://www.plospathogens.org/article/info:doi/10.1371/journal.ppat.1007156#ppat.1007156.t001" target="_blank">Table 1</a>). Asterisk indicates a previously purported miRNA candidate region suggested in the literature [<a href="http://www.plospathogens.org/article/info:doi/10.1371/journal.ppat.1007156#ppat.1007156.ref039" target="_blank">39</a>], which is nonspecific (detected in the negative control JMRV miDGE analysis, lower plot) and can be eliminated by removing low complexity reads (center plot). <b>(B)</b> RNA-seq coverage for the most abundantly mapped HPV in 303 tumors in the TCGA CESC project [<a href="http://www.plospathogens.org/article/info:doi/10.1371/journal.ppat.1007156#ppat.1007156.ref050" target="_blank">50</a>]. Each of the 303 libraries are represented on the X-axis (sorted based on Y-axis value). Y-axis indicates the percentage of the positions in the HPV genome with read mappings. Libraries with > = 50% coverage (213 libraries) were used for subsequent analysis. <b>(C)</b> Percentages of TCGA cervical squamous cell carcinoma (CESC) libraries with miRDeep2 miRNA identifications for each set of reference sequences. Number of libraries examined is 213. <b>(D)</b> Number of unique miRDeep2 miRNA identifications across TCGA CESC libraries for each set of reference sequences. Number of libraries examined is 213. <b>(E)</b> Percentages of Qian <i>et al</i>. [<a href="http://www.plospathogens.org/article/info:doi/10.1371/journal.ppat.1007156#ppat.1007156.ref040" target="_blank">40</a>] libraries with miRDeep2 miRNA identifications for each set of reference sequences. Number of libraries examined is 12. <b>(F)</b> Number of unique miRDeep2 miRNA identifications across Qian <i>et al</i>. libraries for each set of reference sequences. Number of libraries examined is 12. <b>(G)</b> Raw read counts of all small RNAs mapping to the indicated reference sequences for each library from Qian <i>et al</i>.</p