Human foreskin fibroblast (HFF) cells (kindly provided by Dr. Kra-Oz, Haifa, Israel) were cultured in BIO-AMF complete media (Biological Industries, BeitHaemek, Israel) and used to propagate two clinical isolates and the laboratory strain AD169 of HCMV.

HFF cells were infected with two clinical isolates and AD169 laboratory strain, at multiplicity of infection (MOI) of two.

Total RNA was extracted 96 hours post infection using EZ-RNA II isolation kit (Biological Industries, Beit Haemek, Israel), according to the manufacturer's protocol. Aliquots of total RNA were used directly for miRNA quantification using qPCR or subjected to small RNA library construction.

DS was carried out at the Tel-Aviv University Genome High-Throughput Sequencing Laboratory following Illumina’s Small RNA sample preparation protocol v1.5 (TruSeq Small RNA Sample Prep Kits, Illumina) for generating small RNA libraries directly from total RNA and for sequencing. miRNA libraries were sequenced on the Illumina Cluster Station and Genome Analyzer II following the manufacturer’s protocol.

cDNA libraries were prepared using the miScript Reverse Transcription Kit (QIAGEN, Hilden, Germany). Briefly, 1 µg of total RNA (from hCMV infected and uninfected cells) were poly-adenylated by poly (A) polymerase and converted into cDNA by reverse transcriptase with oligo-dT in a single step, at 37^oC for 1 hr. and 95^oC for 5 min using the miScript Reverse Transcription kit (QIAGEN, Hilden, Germany). The cDNA libraries were used in the qPCR analysis for miRNAs using the miScript SYBR Green PCR Kit (QIAGEN, Hilden, Germany) with a forward specific miRNA primer (sequences were taken from the DS analysis, Metabion Martinsried, Germany) and the miScript Universal reverse primer. Forward primers with high G/C content at the 3' end, were extended with one or two 'A' nt at their 3'-end, see Table 1 for list of all primers used in this study. qPCR reactions were performed using the Light Cycler 480 system (Roche Applied Science, Manheim, Germany) and the amplicons were run on a 2% (w/v) agarose gel. Only CMV-miRNAs that were qPCR negative in the uninfected library, and with the correct size of the amplicon were considered as validated.

Sequencing results were first filtered by removing the adaptor sequence and keeping reads with a size of more than 18 nt. The filtered reads were then aligned to the reference Human herpesvirus 5 (HCMV) genome (NC_006273.2) using the BWA program (0.5.8a version downloaded from http://sourceforge.net/projects/bio-bwa/files/), forbidding mismat-ches in the alignment. A genomic map of the virus sequence was built with annotation of its known pre-miRNAs and mature miRNAs (downloaded from miRBase, v.17.0 http://www.mirbase.org/). The results of the BWA were used to assign, to each position in the genome, the reads which start at that position. Based on this, each position was annotated with a read count indicating the number of reads starting at this position (http://www.cs.bgu.ac.il/~vaksler/PROJ/GenomeMapper.htm ).

We inspected the known pre-miRNAs in the genomic map in which only one arm of the precursor contains a known miRNA and checked the number of reads mapped to the other arm. Whenever the read count was above 25, the corresponding sequence was considered as a candidate miRNA and was further checked by qPCR.

We traversed each position in the genomic map in order to identify peaks that are putative start indices of miRNAs. These positions comply with the following criteria: (1) a read count of at least 25, (2) at least 70% of these reads have length ranging between 20-23 nt, and (3) the neighboring positions in the range 6-15 nt on both sides of the considered position all have a read count smaller than 10 (this to guarantee that there are no abundant small RNAs starting in the middle, or too closely above, a putative miRNA, and to take into account the allowed offsets of minus 5 to 5, see Results). In order to further take into account the secondary structure of the potential miRNAs starting at each putative start position, two windows of 100nt in size were considered such that the putative miRNA sequence could participate in either arm of the putative stem. The predicted structures were obtained via MFOLD [17]. The windows whose folding yielded stem-loop hairpins and met the following criteria were considered typical precursors: the sequence with high read count falls in the stem of the hairpin and the free-folding energy is lower than -25 kcal/mol. In addition, the pairing number of the mature sequence within the hairpin should be higher than 14 nts, and the maximal bulge should be less than 5 nts. The putative miRNAs from predicted precursors were further checked by real time PCR (qPCR).

miRNA offset RNAs (moRs) are gaining popularity in the field of small RNAs. We investigated the presence of these 19-20 nt RNAs in our DS data. To find potential moRs we scanned all known and predicted precursors and checked if there were peaks of read counts with starting position approximately 18-20 nt before the 5p-miRNA or at the end of the 3p-miRNA. In addition, we required that the sought moRs be part of a stem in the putative folding.

The aim of our study was to characterize the miRNAome of HCMV. In pursue of known and new miRNAs encoded by HCMV we applied high throughput DS analysis to small RNA prepared from HCMV infected HFF cells.

Among the 16,837,317 reads obtained from the deep-sequencing, 1,015,588 (6%) were perfectly matched to the HCMV genome. These latter reads were assigned to the genomic map of the HCMV using the results from BWA (http://www.cs.bgu.ac.il/~vaksler/PROJ/GenomeMapper.htm). From now on we refer just to these latter reads. 812,477 (80%) of the reads were mapped to the known HCMV miRNAs, 35,727 (4%) of the reads were mapped to new miRNAs, and 7903 (0.7%) of the reads were mapped to moRs. The read count of a specific miRNA includes all the sequences such that their start position deviates from that of the miRBase sequence (our 0-offset) by an offset of -5 to 5 nt. For new miRNAs and moRs, the deviation is from the most abundant sequence (which is our 0-offset).

Fig. (1) illustrates the length distribution of reads that were mapped to miRNAs (black), moRs (dark grey), and other small RNAs (light grey). As known from the literature [18] the length of miRNAs is ~21-22nt. Our DS supports this fact since more than half of the reads mapped to the HCMV miRNAs (known and new) had a length of 22 nt, and ranged from 19 nt to 25 nt. moRs peaked at length of 19-20 nt. The length of the rest of the reads was distributed between 18 to 32nt.

Fig. (1). Length variation of reads that were mapped to miRNAs. Nucleotides length distribution of the percent of reads that were mapped to known and new HCMV miRNAs (in black), moRs (dark grey), and other small RNAs (light grey).

The dominant sequence variants were compared to the consensus data in the miRBase website, and the distribution of the offsets of miRNA start points (on the 5' end) and end points (on the 3' end) from the 0-offset was analyzed. The variation in start points at the 5'-end (Fig. 2A) was smaller than the variation of end points at the 3'-end (Fig. 2B). The frequency of 5' and 3' ends variation was 15% and 65% respectively. The extent of length and offset variations were observed previously in EBV DS analysis as well, and also in human and mouse miRNAs [11,19-21].

Fig. (2). Distributions of 5'-end and 3'-end heterogeneity of HCMV miRNA. A) Distribution of 5'-end heterogeneity of known and new HCMV miRNAs. To analyze the 5'-end heterogeneity, we determined as 0-offset the start position of the miRBAse sequence for known miRNAs. For new miRNAs the 0-offset is the position that had the highest read count. Sequences whose start point was before the 0-offset were assigned a negative offset number, while sequences which started after the 0-offset were assigned a positive offset number. The frequency (% of total reads) of 5' end heterogeneity was calculated by dividing the read count at each offset by the total read count of all offsets at that end. The calculation was based on the 24 mature known and new miRNAs (Tables 2 and 3). B) Distribution of 3'-end heterogeneity of HCMV miRNAs. A similar strategy to A was used to calculate the 3'-end heterogeneity.

HCMV encodes 17 known mature miRNAs from 11 precursors [22]. Sequence information and the read counts for these miRNAs are summarized in Table 2. Interestingly, we found that most of the HCMV miRNAs had many reads, some reaching 100s of thousands (Table 2 and Fig. 3). Thirteen of these miRNAs were found in the sample with a read count of above 3000 reads. The relative abundance of the known miRNAs is depicted in Fig. (3): miR-UL22A-3p had the highest read count (196,688), followed by miR-UL36-5p (116,669) and miR-US25-1-5p (108,624), constituting 23%, 13.7% and 12.7%, respectively, out of all reads mapped to HCMV miRNAs.

Reads for mature miRNAs of mir-UL70-1 and mir-UL148D-1 were not found in our DS. Interestingly, all the reads that were mapped to miR-US4-5p were in +5 offset from the miRBase published sequence.

With respect to already described miRNA biogenesis, we observed that the miRNA/miRNA* designation in the miRBase annotation is consistent with the analysis of our DS data with a single exception of mir-UL22A-1. As opposed to the miRBase and the qPCR (data not shown), we found that mir-UL22A-3p, which has been designated as the star (*) form of this miRNA in the miRBase annotation, had a much higher sequence-read abundance than miR-UL22A-5p (see Table 2).

This alternation in expression of mature miRNAs within their precursor as well as the finding that some miRNAs were not detected by the DS technique while they were detected by another approach, e.g., the qPCR, might be attributed to technical difficulties of the DS such as ligation and preferences in PCR amplification.

Four previously unreported miRNAs, that originated from known pre-miRNAs [3,4], were revealed by the DS: miR-UL112-5p, miR-US4-3p, miR-US5-1-5p, and miR-US5-2-5p, (2865, 1145, 60, and 733 reads, respectively, see Table 3A and Fig. 3). These miRNAs had a relatively low read count. Three of these miRNAs were confirmed by qPCR (Fig. 4), and thus were considered as new HCMV miRNAs. In contrast, miR-US5-1-5p had a different TM (in the TM analysis of the qPCR) and a larger size of amplicons than the expected size of a miRNA.

Fig. (3). HCMV miRNAs read counts and their percentages in the DS. Read counts of all known and new HCMV miRNAs detected by deep-sequencing include all the reads in the offsets minus 5 to 5 of the 0-offset. The percents above the bars (in red) indicate the relative read count of the miRNAs, above threshold of 0.8%.

Fig. (4). qPCR detection of candidate miRNAs encoded by HCMV miRNAs. The known miRNA: miR-UL22A-5p served as positive control. Each experiment was carried out three times, shown are the average Cq and STD for three different experiments on two low passage-number clinical isolates (white and black bars) and one lab-strain (AD169) (gray bars).

Using the method described above, we identified 13 precursors on both strands of the genome, nine of them are known precursors and four are new. The mature miRNAs originating from the latter precursors (see Fig. 5), were further checked by qPCR. Six mature miRNAs from four precursors, mir-US29-1, mir-US22-1, mir-UL69-1 and mir-UL59-1, were confirmed by qPCR (Table 3B and Fig. 4). Furthermore, we also detected these six new miRNAs in RNA extracted from plasma and amniotic fluid samples taken from HCMV infected patients using real-time PCR (data not shown). Fig. (6) depicts the genomic locations of all these new miRNAs, along with the known miRNAs.

Fig. (5). miRNAs and moRs in new (A) and known (B) precursors. The hairpin precursor structures and within them the known (red) miRNAs, new (blue) miRNAs, and moRs (black with frame). The two rightmost columns are the corresponding DS read counts of the 5p and 3p moRs.

Fig. (6). Genomic location of new HCMV miRNAs. Positions of new HCMV miRNAs are shown on fragments of the HCMV genome. The new HCMV miRNAs are in white boxes, and the known miRNAs are in gray boxes.

The correct size of the amplicons of the mature miRNA was analyzed by TM-calling (Light Cycler 480) and on agarose gel.

mir-UL59-1 and mir-UL69-1 are encoded within an intergenic region between genes UL57 and UL69. Even though the 3p arm of mir-UL59 was detectable by the DS, we could not confirm it by qPCR. Both the 5p arms of these two miRNA precursors are readily detectable in HFF cells infected by clinical isolates and strain AD169 (Fig. 4) as well as in plasma and amniotic fluid specimens obtained from HCMV infected patients (data not shown).

The mir-US22-1 precursor, which contains two mature miRNAs, is transcribed from the complementary strand to gene US22, and thus the 5p and 3p mature miRNAs are likely candidates to down regulate US22 expression. (This phenomenon is similar to that observed in other previously known miRNAs (UL70-1, UL112-1, UL148D, US33-1) [3, 4].

mir-US29 that has been predicted before [3], is located within the US29 gene and is transcribed antisense to the previously known mir-US33-1 (Fig. 6). The phenomenon that two precursors are transcribed from both strands of a genomic region has been previously reported in Drosophila [23], HSV-1 [14] and in KSHV [10].

In order to further strengthen the validity of the new HCMV miRNAs that were found by DS and confirmed by qPCR, we looked for sequence homology between the pre-miRNAs of HCMV and the Chimpanzee cytomegalovirus (CCMV), the closest relative of HCMV in herpesvirus β family [24]. We could find similarity of 100% and 69%, respectively, between pre-miRNAs mir-US29-1 and mir-US22-1 of HCMV and the CCMV sequences. Both corresponding sequences in CCMV could be folded into hairpin secondary structures (data not shown). No homologous sequence for mir-UL59 and mir-UL69 were found in the CCMV genome that could fold into a hairpin secondary structure.

We also checked sequence homology between the mature miRNAs of HCMV and sequences of CCMV. Both miRNAs: miR-US29-5p and miR-US29-3p were fully conserved in CCMV. The seed of miR-US5-2-5p is fully conserved in CCMV, while the seeds of miR-US4-3p, miR-US22-5p, miR-US22-3p and miR-UL112-5p show a single nucleotide mismatch with CCMV.

Up to this date no moRs of HCMV have yet been reported. In this work we identify 11 moRs, denoted by moR-UL22A-1-5p, moR-UL22A-1-3p, moR-UL112-5p, moR-US4-1-5p, moR-US4-1-3p, moR-US5-1-5p, moR-US5-2-5p, moR-US25-1-5p, moR-US33-1-5p, moR-US22-1-5p, and moR-US29-1-5p. In Fig. (5A, B) we present the moRs of HCMV within the secondary structure of their precursors. We also list in the figure the number of reads of the moRs. The number of reads of moRs ranged between tens to thousands. The high abundance of some of the moRs suggests that they might have important biological functions. In order to further investigate these small RNAs we selected four of the most abundant moRs: moR-US-5-2-5p, moR-UL22A-1-5p, moR-UL112-5p, and moR-US4-5p for qPCR expression analysis in our two low-passage clinical isolates and strain AD169. The qPCR detected three out of the four moRs: moR-US-5-2-5p, moR-UL112-5p, and moR-US4-5p (Fig. 7) in all the three HCMV isolates. Size analysis of qPCR products by TM calling and gel electrophoresis confirmed their size. In spite of its high read count, the qPCR failed to detect moR-UL22A-1-5p. The fact that HCMV moRs are expressed in authentic infection, and their consistent expression in various isolates, strongly suggests that these small RNAs have important biological roles in HCMV infection.

Fig. (7). qPCR detection of HCMV three moRs in three HCMV strains. Shown are: moR-US4-5p (black bars), moR-UL112-5p (white bars), and moR-US5-2-5p (gray bars). Each experiment was carried out three times, shown are the average Cq and STD between experiments.