Next‐generation sequencing of miRNAs in clinical samples of Epstein–Barr virus‐associated B‐cell lymphomas

Abstract Epstein–Barr virus (EBV) encodes 49 microRNAs (miRNAs) in the BART and BHRF1 regions of its genome. Although expression profiles of EBV‐encoded miRNAs have been reported for EBV‐positive cell lines and nasopharyngeal carcinoma, to date there is little information about total miRNA expression, including cellular and viral miRNAs, in the primary tumors of EBV‐associated B‐lymphoproliferative disorders. In this study, next‐generation sequencing and quantitative real‐time reverse transcription‐PCR were used to determine the expression profiles of miRNAs in EBV‐infected cell lines and EBV‐associated B‐cell lymphomas, including AIDS‐related diffuse large B‐cell lymphoma (DLBCL), pyothorax‐associated lymphoma, methotrexate‐associated lymphoproliferative disorder, EBV‐positive DLBCL of the elderly, and Hodgkin lymphoma. Next‐generation sequencing revealed that EBV‐encoded miRNAs accounted for up to 34% of total annotated miRNAs in these cases. Expression of three miR‐BHRF1s was significantly higher in AIDS‐related DLBCL and pyothorax‐associated lymphoma compared with methotrexate‐associated lymphoproliferative disorder and EBV‐positive DLBCL of the elderly, suggesting the association of miR‐BHRF1s expression with latency III EBV infection. Heat map/clustering analysis of expression of all miRNAs, including cellular and EBV miRNAs, by next‐generation sequencing demonstrated that each EBV tumor, except methotrexate‐associated lymphoproliferative disorder, formed an isolated cluster. Principal component analysis based on the EBV‐encoded miRNA expression showed that each EBV tumor formed a distinguished cluster, but AIDS‐related DLBCL and pyothorax‐associated lymphoma formed larger clusters than other tumors. These data suggest that expression of miRNAs, including EBV‐encoded miRNAs, is associated with the tumor type and status of virus infection in these tumors.


Introduction
Epstein-Barr virus (EBV) is associated with the pathogenesis of malignant tumors, including lymphoma, gastric cancer, and nasopharyngeal carcinoma [1,2]. These EBVassociated tumors have different patterns of EBV-encoded mRNA expression, which are categorized as latency types I, II, and III. Type III latency is characterized by the expression of all EBV latent genes, including EBV-encoded nuclear antigens (EBNAs) and latent membrane proteins (LMPs). In EBV-associated malignancies with type I and II latency, the expressions of EBV genes are restricted by the activity of a specific promoter. While EBNA1 and EBV-encoded small RNAs (EBERs) are expressed in almost all types of EBV-infected cells, the expressions of LMPs and other EBNAs are downregulated in EBV-positive Burkitt lymphoma and gastric cancer. Different expression profiles are thought to contribute to the pathogenesis of EBV-associated malignancies [2]. LMP1 and EBNA2 are expressed in malignancies of type III latency, such as opportunistic EBV-associated lymphoma.
In addition to viral miRNAs, cellular miRNAs also play important roles in the oncogenesis of B-cell lymphomas [22]. miR155 is accumulated in B-cell lymphoma cells and plays a crucial role in the growth of B-cell lymphoma [23,24]. miR155 has been shown to suppress activation-induced cytidine deaminase-mediated MYC-IGH translocation in Burkitt lymphoma [25]. miR19a and miR19b in Burkitt lymphoma cells are a direct transcriptional target of C-MYC [26]. EBV infection also induces expression of oncomiRs, such as miR-21 and miR-146a, in B-cells [27][28][29].
To clarify the roles of cellular and EBV-encoded mi-RNAs in vivo, it is important to determine the expression profile of total miRNAs in primary tumors of EBVassociated lymphoma. To date, the expression profiles of EBV-encoded miRNAs have been investigated in primary NPC, nasal NK/T-cell lymphoma, gastric cancer samples, and LCLs using real-time RT-PCR [10,11,14,[30][31][32][33][34][35][36][37], but there is little information about EBV-encoded miRNA expression in primary tumors of EBV-associated B-cell lymphoma. Next-generation sequencing (NGS) is a powerful tool for investigating the expression of small RNAs, including cellular and viral miRNAs. The expression profile of cellular and viral miRNAs has been reported in EBVpositive diffuse large B-cell lymphoma (DLBCL) of the elderly and Burkitt lymphoma cases by NGS, DNA microarray, and real-time PCR analysis [33,[37][38][39][40]. However, no report has compared the miRNA expression profiles among different histological categories of EBV-associated B-cell lymphomas. In this study, the expression of cellular and EBV-encoded miRNAs was investigated using NGS and quantitative real-time RT-PCR analysis in EBV-positive cell lines and clinical samples of EBV-associated B-cell lymphomas.

Materials and Methods
Samples Studies using human tissue were performed with the approval of the Institutional Review Board of the National Institute of Infectious Diseases (Approval No. 271 and 272) and Tokyo Metropolitan Komagome Hospital (Approval No. 628). All clinical samples were collected from Japanese patients in Japan. All data in this study were analyzed anonymously. Seventeen frozen samples of EBV-associated diseases were examined by NGS. Samples included three AIDS-related DLBCL (ARL), four pyothorax-associated lymphoma (PAL), five methotrexaterelated lymphoma (MTX), three EBV-positive DLBCL of the elderly (ELD), and two Hodgkin lymphoma (HL) cases. In addition, 18 ARL, four PAL, four MTX, five ELD, and six HL cases were analyzed by real-time RT-PCR. The four ARL case samples that were analyzed by real-time RT-PCR were formalin-fixed paraffin-embedded (FFPE) tissues. Histological diagnosis of lymphoma was based on the 4th edition of World Health Organization classification of lymphoma [41]. All cases were confirmed by immunohistochemistry or RT-PCR, suggesting latency III. Cases of MTX and ELD were negative for LMP1, but some of them were not examined. HL cases were positive for LMP1 and negative for EBNA2 in Reed-Sternberg cells, suggesting latency II. Histological subtypes of DLBCL included centroblastic and immunoblastic variants, but some cases were not distinguishable histologically. Six EBV-negative DLBCL tumors from immunocompetent patients were also examined. In addition to the tumors from clinical cases, one EBV-transformed EBV LCL (spontaneous transformation), four EBV (strain B95-8)transformed LCLs, and two EBV-negative and KSHV-positive primary effusion lymphoma cell lines, TY-1 and BCBL-1, were examined.

RNA extraction
Total RNA, including miRNA, was extracted from frozen or paraffin-embedded samples and cell lines using the High Pure miRNA purification kit (Roche Molecular Biochemicals, Indianapolis, IN) and subsequently treated with Turbo DNase (Ambion, Austin, TX), all according to instructions from the manufacturers.

NGS
Small RNA libraries were established with the TruSeq small RNA kit (Illumina, San Diego, CA) from 18-to 35-nucleotide cDNAs using 5 μg of DNase-treated total RNA. Small RNA was sequenced using the Illumina Genome Analyzer IIx platform with cBot-SR v2 and SBS 36 cycle v5 kit. Sequence reads were analyzed with CLC Genomics Workbench (version 7.5.1, CLC bio, Aarhus, Demark). miRBase release 21 was used as the miRNA database (http://www.mirbase.org/).

Real-time RT-PCR for miRNA
Total RNA (100 ng) was reverse-transcribed using the miScript Reverse Transcription Kit from Qiagen (Valencia, CA). Real-time reverse transcription polymerase chain reaction (RT-PCR) for the quantification of 42 EBVencoded miRNAs and human cellular miRNAs miR16 and miR21 was carried out with the miScript PCR system (Qiagen). Ratios of the copy numbers of EBV-encoded miRNA to miR16 were calculated as follows: ratio of target miRNA to miR16 = 2 Ct of miR16 /2 Ct of target (Ct = cycle threshold). Each Ct was determined by results of EBVnegative controls for EBV-encoded miRNA and nontemplate controls for miR16 in each experiment. In addition, stem-loop real-time RT-PCR (Taqman microRNA assay, Applied Biosystems, Foster City, CA) was also performed for comparison with the miScript PCR system. Synthesized miR21 with deletion or addition of single or double nucleotides in the 3′ end (Fasmac, Tokyo, Japan) was examined.

Heat map and cluster analysis
Heat map images were produced using TreeView and Cluster software by Michael Eisen, University of California at Berkeley (http://rana.lbl.gov/EisenSoftware.htm) [42].

Principal component analysis
Principal component analysis (PCA) was performed on normalized data of ratios of EBV-encoded miRNA copy number to miR16 copy number using the PCA function of SPSS software (IBM, Armonk, NY). The first three principal components were used to produce twodimensional and three-dimensional plots. Resolution and convex hulls in three-dimensional plots were calculated and plotted using MATLAB (MathWorks, Natick, MA), as reported previously [31]. Volumes of convex hulls were measured by MATLAB ConvexHull function after Delaunay-triangulation.

Accession number
The annotated miRNAs detected in cell lines and clinical samples by NGS in this study are registered as accession number DRA002823 in the DDBJ Sequence Read Archive.

Statistical analysis
Mann-Whitney U-test was used for nonparametric twogroup comparison (Graph Pad Prism 5, GraphPad Software, La Jolla, CA).

Small RNA libraries
Small RNA libraries were established from an EBVtransformed LCL and 17 clinical samples of EBV-associated lymphoproliferative disorders to determine the miRNA profiles. Although more than 100,000 reads were sequenced in each sample by NGS, unannotated rates ranged from 13.91% to 65.61% (Table S1). Noncoding RNA reads ranged from 14.43% to 79.38% and miRNA reads ranged from 1.16% to 71.66%. A total of 688 cellular and EBV miRNAs were annotated from all small RNA libraries (Table S2). EBV-encoded miRNAs ranged from 0.03% to NGS of miRNA in EBV B-cell Lymphoma K. Sakamoto et al. 34.7% of all annotated miRNA reads, and their percentages varied among cases (Table S1).
Since all reads matching to pre-miRNAs were counted as miRNAs in this system, the read counts contained both 5p and 3p mature miRNAs but also immature miRNAs. We also investigated the profile of coverage in each pre-miRNA. Coverage profiles in pre-miRNAs demonstrated frequent deletion or addition of a single nucleotide in the 3′ end of mature miRNAs in both cellular and EBV-encoded miRNAs in a spontaneous transformed LCL (Fig. 4). Among 20 miRNAs picked up randomly from cellular and viral miRNAs, the frequencies of changes in the 3′ end of mature miRNAs were not different between     those of cellular and viral miRNAs (data not shown). In addition, frequencies of changes in the 3′ end of mature miRNAs were not different between miRNAs frequently expressed in clinical samples of EBV-associated lymphomas and those in LCL (data not shown).

Expression profile of EBV-encoded miRNAs using real-time RT-PCR
To confirm the expression profile of EBV-encoded mi-RNAs, we performed real-time RT-PCR of these EBVencoded miRNAs using additional cases of EBV-associated lymphoma (18 ARL, 4 PAL, 4 MTX, 5 ELD, and 6 HL cases). Since NGS demonstrated frequent deletion or addition of a single nucleotide in the 3′ end of mature miRNA, we first compared the miScript PCR system and stem-loop real-time RT-PCR by examining synthesized miR21 that contained a deletion or addition of one or two nucleotides in the 3′ end. The miScript PCR system succeeded in detecting miR21 with deletions and additions, whereas stem-loop real-time PCR detected miR21 with deletion or addition with less than one-tenth sensitivity (Fig. 5A). Thus, the miScript PCR system was employed to detect 5p and 3p mature miRNAs separately. Real-time RT-PCR analysis revealed a similar profile of EBV miRNAs between FFPE tissues and frozen samples   (Fig. 5B). In addition, cycle thresholds of miR16 were proportional to those of miR21 in 12 representative clinical samples (Fig. 5C); therefore, the copy numbers of EBV-encoded miRNAs can be normalized by those of miR16 in each sample. A total of 1 × 10 5 copies of synthesized miR21, miR21 with deletion (miR21-1), miR21 with addition of single nucleotide (miR21 + 1), and miR21 with addition of two nucleotides (miR21 + 2) in the 3′ end of mature miRNA were examined with the miScript PCR system (miScript SYBR green real-time PCR) and stemloop real-time RT-PCR (Taqman microRNA assay). Error bars indicate standard error for triplicate analyses. (B) Comparison of formalin-fixed paraffinembedded (FFPE) and frozen samples by real-time RT-PCR. miR-BARTs and miR-BHRF1s were quantified by the miScript PCR system. FFPE and frozen samples were obtained from xenotropically inoculated lymphoma tissues in severe combined immunodeficiency mice. (C) Cycle thresholds (Ct) for miR16 and miR21. The copy numbers of two cellular miRNAs, miR16 and miR21, were measured and plotted in 12 representative clinical samples using the miScript PCR system. Linear approximation line (broken line) and correlation coefficient (R 2 ) are indicated.   Real-time RT-PCR revealed that the pattern of EBVencoded miRNA expression in ARL is similar to PAL (Fig. 6). miR-BHRF1-1 showed higher copy numbers in ARL and PAL than in MTX, ELD, and HL (Mann-Whitney U-test, P < 0.01). In HL, we detected miR-BARTs 13-3p and 19-3p, with no or low level expression of the other miRNAs. In addition, total expressions of EBV-encoded miRNAs were low in MTX and HL cases compared with other cases (Mann-Whitney U-test, P < 0.01).

Heat map, cluster analysis, and PCA
To determine if EBV-associated lymphoma cases were categorized based on expression patterns of miRNAs, we performed cluster analysis on our data set using heat maps and PCA. Heat map and cluster analysis based on the data by NGS demonstrated that ARL, PAL, HL, and ELD formed an individual cluster based on the data of total miRNA expression, including both cellular and viral miR-NAs (Fig. 7). MTX cases belonged to various clusters and did not form a cluster. PCA based on the data by EBVencoded miRNA expression by real-time RT-PCR analysis showed small volumes of cluster in three-dimensional plots of PC1-3 corresponding to LCL (cluster volume = 0.023), HL (1.62 × 10 −5 ), MTX (8.63 × 10 −4 ), and ELD (0.029) samples, whereas ARL (1.426) and PAL (2.030) formed larger clusters (Fig. 8). PCA showed that each EBV tumor, other than ARL and PAL, formed a distinguished cluster, but samples of ARL and PAL distributed broadly in the three-dimensional plots of PC1-3 (Fig. 8D).

Discussion
In this study, we describe the expression profiles of cellular and viral miRNAs in primary tumors of EBVassociated lymphoma. High expression of miR21 is common among many types of tumors, including B-cell lymphoma [46][47][48]. In addition, miR21 was identified as a useful marker for primary diffuse large B-cell lymphoma of the central nervous system in the cerebrospinal fluid [49]. miR143 has a possible role in angiogenesis in tumors, suggesting an association between miR143 expression and lymphoma growth [43]. Downregulation of miR143 was demonstrated in B-cell lymphoma cell lines [50]. In this study, NGS revealed 0.008% of total annotated reads of miRNA in an EBV-transformed LCL (Fig. 2). Low expression of miR143 in EBV-positive cell lines was previously reported, such as 0.003% in Jijoye cells, 0.0004% in LCL35, 0.0002% in LCL-BAC, and 0.0007% in SDLCL in total annotated reads of miRNAs [7,51]. On the other hand, this study demonstrated that miR-143 comprised 6.07% of total miRNAs on average in clinical samples. This observation indicates a difference in miRNA expression between cell lines and clinical primary cases. We also confirmed high copy number of miR155 in EBV lymphoma in this study. Lower copy number of miR155 in ELD and HL than ARL and PAL may be associated with the presence of nontumor cells and inflammatory cells in the samples (Fig. 2). In particular, tumor cells are observed at low frequency in the lymph node with HL, whereas almost all cells were tumor cells in tumor tissues of ARL and PAL (Fig. 1).
Several reports have described expression levels of EBVencoded miRNAs using stem-loop RT-PCR [11,52,53] and NGS [33] in EBV-associated diseases. In our study, NGS identified frequent variants of mature miRNAs with 3′ deletion and addition (Fig. 4), which the stem-loop Figure 7. Heat map and cluster analysis. Heat map and cluster analysis of miRNA expression using data of NGS. Individual patient samples are shown in columns and miRNAs including both cellular and EBV miRNAs in rows. Expression measured by NGS is displayed in red color and green color depending on expression above or below median expression level.

Low expression
High expression NGS of miRNA in EBV B-cell Lymphoma K. Sakamoto et al.
RT-PCR was unable to detect (Fig. 5). This is important information for detecting miRNAs using the stem-loop primer, because the stem-loop RT-PCR may fail to detect large portion of miRNAs with such a deletion or addition of nucleotide in the 3′ region.
In this study, the miScript PCR system failed to detect almost all miRNAs encoded by BART cluster 2 in the LCL transformed by B95-8 EBV (Fig. 6). These results were reasonable, as cluster 2 is missing in B95-8, which has been widely used in the establishment of LCLs. miR-BARTs in cluster 2 were not detected in LCLs transformed by B95-8 EBV, but were detected at high levels in spontaneously transformed LCLs, ARL, and other clinical samples. A recent study revealed that at least 8 miRNAs in BART cluster 2 (miR-BARTs 15, 18-3p, 7-5p, 10-3p, 10-5p, 11-3p, 13-5p, and 14-5p) were identified as associated with latency type III infection of EBV [31]. Our data confirmed this, in part. For example, miR-BARTs 15, 10-3p, 11-3p, and 14-3p were expressed at higher levels in ARL than other tumors. However, expression levels of miR-BARTs 7-5p, 18-5p, and 13-5p were low in all samples. These observations are consistent with previous studies investigating lymphoma samples [31,33], although others reported no impact of miR-BARTs on LCL growth or survival in vitro [15]. As reported previously [7,51,54], miR-BHRF1 expression is associated with latency III infection of EBV. Our results by NGS and real-time RT-PCR demonstrate relatively high copies of miR-BHRF1-1, -2, and -3 expressions in ARL and PAL, both of which are latency III infection of EBV. Thus, our data suggest that miRNA expression differs between tumors and may be associated with EBV latency, in part.
PCA showed a broad distribution of ARL and PAL samples, but close clusters were formed corresponding to HL, MTX, ELD in two-and three-dimensional plots (Fig. 8), suggesting that EBV-encoded miRNA expression was deregulated in ARL and PAL. One study reported that EBVencoded miRNAs were strictly regulated in normal EBV-positive nontumor B-cells, and that deregulation of EBV-encoded miRNA expression was observed only in  [31]. The results of heat map and PCA in this study suggest that EBV-encoded miRNA expression was partially regulated in HL, MTX, and ELD but drastically deregulated in ARL and PAL. The most probable difference in clinical manifestations between ARL/PAL and others is the immune status of the host. While the mechanism of viral miRNA expression is unknown, the immune status of the host may affect viral miRNA expression. Further studies are required to reveal details of the association between viral miRNA expression and host immunity.