EBV miRNAs are potent effectors of tumor cell transcriptome remodeling in promoting immune escape

The Epstein Barr virus (EBV) contributes to the tumor phenotype through a limited set of primarily non-coding viral RNAs, including 31 mature miRNAs. Here we investigated the impact of EBV miRNAs on remodeling the tumor cell transcriptome. Strikingly, EBV miRNAs displayed exceptionally abundant expression in primary EBV-associated Burkitt’s Lymphomas (BLs) and Gastric Carcinomas (GCs). To investigate viral miRNA targeting, we used the high-resolution approach, CLASH in GC and BL cell models. Affinity constant calculations of targeting efficacies for CLASH hits showed that viral miRNAs bind their targets more effectively than their host counterparts, as did Kaposi’s sarcoma-associated herpesvirus (KSHV) and murine gammaherpesvirus 68 (MHV68) miRNAs. Using public BL and GC RNA-seq datasets, we found that high EBV miRNA targeting efficacies translates to enhanced reduction of target expression. Pathway analysis of high efficacy EBV miRNA targets showed enrichment for innate and adaptive immune responses. Inhibition of the immune response by EBV miRNAs was functionally validated in vivo through the finding of inverse correlations between EBV miRNAs and immune cell infiltration and T-cell diversity in BL and GC datasets. Together, this study demonstrates that EBV miRNAs are potent effectors of the tumor transcriptome that play a role in suppressing host immune response.


Introduction
The Epstein Barr Virus (EBV) is a ubiquitous gammaherpesvirus (γHV) that establishes lifelong infections in over 90% of the world's population. Since its discovery as the major etiological agent of endemic Burkitt's lymphoma (BL) [1], EBV has been causally associated with other malignancies including NK/T cell lymphoma [2], diffuse large B cell lymphoma [3], Hodgkin's lymphoma [4], nasopharyngeal carcinoma [5] and gastric carcinoma (GC) [6]. Arising from infected founder cells, these tumors maintain a dependence on the virus as they progress [7].
γHVs utilize two distinct strategies, referred to as lytic and latent replication, to expand the infected host cell population [8]. The lytic replication program is a process utilized by all viruses to replicate and package their genetic content for spread through cell-to-cell transfer. This strategy produces a large number of virions but it destroys the host cell and triggers a strong local immune response [9]. The viral latency program, the phase most closely linked with cancer, is associated with the expression of a small number of genes that support the growth and health of the infected cell [10][11][12][13][14][15][16] and by extension, the growth and health of its viral occupant. In this phase of the virus infection cycle, viral and host genomes are replicated concordantly and distributed to daughter cells [17], resulting in an expansion of the infected cell population that is independent of virion production. This form of intracellular replication minimizes immune reactivity, allowing the virus to discretely double the pool of infected cells at each mitotic cycle.
In EBV-associated BLs and GCs, the latency gene expression program is especially restrictive, with the only ubiquitous and consistently expressed viral protein being the replicative factor, EBNA1 [18]. More abundantly expressed are a group of viral noncoding RNAs that allow the virus to modulate the host cell environment without eliciting a strong adaptive immune response. These include two short noncoding RNAs (EBER1 and EBER2) [19], a long noncoding RNA (RPMS1) [20], circular RNAs [21,22], and the viral miRNAs [23].
MiRNAs function by targeting complementary mRNAs for destruction [24]. By interacting with an average of 90 unique transcripts[25], a single miRNA can have a substantial impact on the cellular transcriptome [26][27][28][29][30]. The genomes of EBV, KSHV, and MHV68 potentially encode 44, 25, and 28 miRNAs [31][32][33][34], theoretically endowing these viruses with the capacity to control nearly every pathway in the cell. Several previous studies have used broad-scale approaches such as AGO-CLIP to uncover viral miRNA targets [35][36][37], and some of these miRNA-target interactions have since been shown to interrupt apoptosis [38], prevent reactivation [39], and block interferon signaling [40]. The AGO-CLIP approach, however, requires bioinformatic inferences to determine miRNA-mRNA pairings, precluding analysis of binding efficacy and limiting detection to canonical interactions. Here we used a modified version of a more comprehensive and quantitative approach, Crosslinking, Ligation, and Sequencing of Hybrids (CLASH) [41], called qCLASH [42], to broadly uncover bona fide targets of EBV miRNAs. Through integration of mRNA, miRNA, and hybrid abundances, we examine the global binding properties and targeting efficacies inherent to γHV miRNAs. Assessing pathway enrichment of the highest targeting efficacy EBV miRNA-target interactions shows enrichment for innate and adaptive immune responses. Utilizing clinical BL and GC datasets, we explore the roles of EBV miRNAs in primary tumors and demonstrate a role for EBV miR-NAs in dampening the immune response to viral infection and oncogenesis.
Aligned miRNA sequencing reads were converted to FASTQ format using the following command: samtools fastq -F 0x900 $bamfile Raw miRNA sequencing reads were aligned to an index of mature human and EBV miRNA sequences (miRbase v22) via bowtie with the following command: bowtie $mir_index -n 3 -m 10 --best --strata -p 10 -S $mir_fastq UV Crosslinking qCLASH was performed as previously described 1 , with 3 biological replicates processed for each cell line. SNU719 cells were trypsinized, resuspended in RPMI 1640 plus 10% FBS, then washed twice with PBS. Akata cells were washed twice with PBS. Both cell lines were resuspended in 10 mL PBS then crosslinked with 250 nm λ UV light for a total exposure of 600J/ cm 2 . 50 million cells were transferred to a 1.5 mL tube, pelleted at 800 x g for 5 minutes, and stored at -80˚C before processing as described below.

Cell lysis
Cell pellets were thawed and resuspended in 500 μl cell lysis buffer (

Elution and RNA extraction
Bead complexes were incubated in 100 μl of elution buffer (100 mM NaHCO 3 , 1% SDS) for 15 minutes at 25˚C on a 1,400 rpm shaker. After spinning at 1,000g for 1 minute, the elution buffer was transferred to a new tube. An additional 100 μl of elution buffer was added to the bead complexes, elution was repeated and the eluates were then combined. To improve RNA phase separation of crosslinked RNA-AGO complexes, proteins were digested using proteinase K. 10 μl of proteinase K (Roche, catalog no. 03115887001) in 40 μl of proteinase K buffer (100 mM Tris-HCl (pH 7.5), 50 mM NaCl, 10 mM EDTA) was added to the 200 μl of eluate and samples were incubated at 37˚C for 20 min. Samples were then subjected to phenol-chloroform extraction to purify the RNA.

Targeting efficacy
Targeting efficacy for each miRNA-mRNA hybrid was calculated using the following equation:

AGO RNA immunoprecipitation
Approximately 5x10 6 SNU719 cells per reaction were lysed for 15 minutes on ice in 1 mL NP40 buffer (150 mM NaCl, 50mM Tris-HCl pH 8.0, 1% NP40 (Millipore Sigma, catalog no. 56741)) with 20uL of RNaseOUT (ThermoFisher Scientifc, catalog no. 10777-019) and 10 uL of Protease Inhibitor Cocktail (Cell Signaling, catalog no. 5871). Lysates were spun at 20,000g for 15 minutes at 4C. Supernatants were transferred to a new tube. For each reaction, 30uL of Protein G magnetic beads (New England Biotechnology, catalog no. S1430S) were washed 3 times in 0.1M sodium phosphate and split into two tubes (one for preclearing lysates, the other for the pulldown of antibody-antigen complexes). For preclearing, the equivalent of 15uL of beads were added to lysates and rotated at 4C for 30 minutes. Tubes were put on a magnetic column and supernatants were transferred to a new tube. Either Mouse anti-pan Ago clone 2A8 (1:100 dilution; EMD Millipore, catalog no MABE56) or mouse anti-HA (1:50 dilution; Santa Cruz Biotechnology, catalog no. sc-7392) antibodies were added to precleared lysates and the tubes were rotated for 2h at 4C. The remaining Protein G beads were added, and tubes were rotated for 30m at 4C. Tubes were placed on a magnetic rack and supernatants were discarded. Beads were washed 3x with NP40 buffer then reconstituted in 150 uL of TE buffer (10mM Tris-HCl pH 8.0, 1 mM EDTA). 3 uL of proteinase K was added to each tube and tubes were incubated at 55C for 30 minutes. RNA was isolated via Phenol-Chloroform extraction followed by ethanol precipitation. microRNA RT-qPCR microRNA quantification was performed using a two-tailed RT-PCR approach [49]. RNA was reverse transcribed using M-MLV reverse transcriptase (Millipore Sigma, catalog no. M1302) according to the manufacturer's protocol, using pooled microRNA-specific RT primers at a final concentration of 1nM (see S8 Table for primer sets). qPCR using 1 μl of the resulting cDNA was performed using iQ SYBR green Supermix (Bio-Rad, catalog no. 170-8882) on a Bio-Rad CFX96 instrument using the primers listed in S8 Table. Amplification began with an initial denaturation step at 95˚C for 3 minutes, followed by 40 cycles of 15 second denaturation at 95˚C and 1 minute of annealing and extension at 53-60˚C, depending on the predicted melting temperature of each primer set.

Pathway analysis
CLASH derived mRNA targets for each virus and corresponding host were ranked by targeting efficacy in each cell line. Protein-protein interaction (PPI) and expression correlation networks were assembled for each of the top 20 targets using STRING[53] with default parameters, limited to 20 direct and 10 additional node proteins per queried "seed" gene. The resulting genes were explored for pathway enrichment with the enrichR[54] API, interrogating the Reactome library of pathways. Pathway enrichment was considered significant if the adjusted P-value < 0.05. All significant pathways were assigned to each seed gene.

High expression of viral miRNAs in EBV-associated tumors
BLs and GCs are aggressive tumors with distinct etiologies. Pediatric BL is a B-cell malignancy endemic to sub-Saharan Africa that is characterized by the t(8;14)(q24;q32) MYC:IGH translocation [55,56]. GCs encompass a diverse group of epithelial tumors [57] originating from the stomach lining that have a broader epidemiology [58]. Despite overt differences in pathology and etiology, nearly all endemic BLs (over 90%)[47] and a subset of GCs (~10%) [6] are causally infected with EBV. As a first assessment of the viral contributions to these tumors, we performed principal component analysis (PCA) of 86 BL[47] (70 EBV+) and 235 GC [59] (24 EBV +) cell transcriptomes. Applying logistic regression to each of the top 10 (BL) or 15 (GC) principal components revealed that EBV status is a major distinguishing clinical covariate in both tumor types (Figs 1A and S1), indicating that EBV is likely a key determinant in shaping the tumor transcriptome.
To investigate underlying features of EBV that contribute to the tumor phenotype, we first carried out a quantitative assessment of the viral transcriptomes in EBV-associated BLs and GCs to identify candidate viral effectors in the natural tumor setting. Unlike lytic conditions, where 20% of expressed poly-adenylated RNAs in the cell are of viral origin (S2 Fig), the sum of expressed latency transcripts in BLs and GCs does not exceed 0.02% of cellular RNAs (Figs 1B and S3), with the majority of these being non-coding. Further, some of the EBV transcript signal in this analysis is derived from lytic gene expression occurring in a minor population of tumor cells through occasional sporadic reactivation. Therefore, long RNA latency viral gene expression represents only a minor portion of the overall transcriptome in these datasets. We were unable to assess expression of the highly abundant non-polyadenylated small non-coding viral EBER1 and EBER 2 transcripts because the GC dataset was polyA-selected and because they are partially lost following RNA size selection in the ribodepletion-derived BL dataset (which results in unreliable quantifications across samples).
Housed within the most abundantly expressed viral polyadenylated RNA, the long noncoding RNA (lncRNA), RPMS1, are 20 densely clustered intronic pre-miRNAs encoding at least 31 mature miRNAs (Figs 1C, S4, and S5). In contrast to the limited contribution of viral long RNAs to the tumor transcriptomes, the cumulative expression of viral miRNAs is remarkably high, exceeding host totals by as much as 3-fold ( Fig 1D). This is particularly relevant because miRNA function is critically dependent on sufficient miRNA expression to saturate a high fraction of each target transcript. The importance of miRNA expression levels on function is supported by previous studies that showed that poorly expressed miRNAs have little discernable biological activity [60] whereas sufficient expression of an individual miRNA can have a marked impact on tumor biology [61]. The overwhelming expression of 31 viral miR-NAs ( Fig 1D and S1 Table) in EBV-positive BLs and GCs supports their likely relevance in modulating the tumor transcriptome landscape.  [59] (GC). Long fraction sequences were aligned to the combined human (GENCODE GRCh38.p13) [44] and EBV [45] transcriptomes. Short fraction sequences were aligned to human and viral miRNA sequences (miRbase v22

EBV miRNAs are over-represented in mRNA-bound RNA-induced silencing complexes (RISCs)
With the dependence on miRNA loading into RISC for target destabilization, the level of miR-NA-RISC association is a predictor of miRNA targeting [62]. To assess RISC association characteristics of EBV and host miRNAs, we performed a modified version [42] of Crosslinking, Ligation, And Sequencing of Hybrids (CLASH) [63], referred to as qCLASH [42], in EBV+ cell lines modeling BL (Akata) and GC (SNU719) (S6A Fig). In CLASH, each AGO-bound miRNA-mRNA pair is ligated and then sequenced as a contiguous read, reproducibly resolving both the diversity of miRNA targets as well as the relative abundance of each interaction (S6B Fig).
In conjunction with our qCLASH analyses, we also performed standard small RNA fraction sequencing to assess the underlying expression of each miRNA in these cell lines. Notably, viral miRNA expression in Akata and SNU719 cell lines was substantially lower than the levels observed in primary tumors (Figs 2A and 2B, left, and S7). Nevertheless, this finding is in-line with previous studies showing that EBV miRNA levels are low in cell lines (including SNU719 cells) but increase markedly upon passage in immunocompromised mice [64,65], likely due to an enhanced reliance on viral miRNAs in the in vivo setting. Despite the modest expression of viral miRNAs in in vitro cultured SNU719 and Akata cells, we found a remarkably high representation of EBV miRNA-mRNA hybrids (Fig 2A and 2B, left). This observed enrichment was not attributable to a limited number of individual viral miRNAs with unusually high binding efficiencies, but instead was a characteristic that was broadly attributable to the bulk of EBV miRNAs (Fig 2A and 2B, right). These results indicate that in addition to exhibiting high in vivo expression, EBV miRNAs evolved with feature(s) that enhance RISC formation, supporting the possibility of distinctly productive viral miRNA targeting.
To determine whether enhanced representation of viral miRNAs in RISC is a conserved feature of γHVs, we supplemented data from our previously published KSHV [42] and MHV68 [66] qCLASH studies with small RNA sequencing to similarly compare KSHV and MHV68 hybrid levels to their respective miRNA expression. Similar to our findings with in vitro EBV miRNA expression, KSHV miRNAs were found to be expressed at low levels in KSHV infected LTC-TIVE cells (S8A Fig). Like our observations with EBV miRNAs, KSHV miRNAs exhibited marked enrichment in RISCs that was distributed across a spectrum of KSHV miRNAs (S8A Fig). In contrast, an aggregate assessment of MHV68 miRNAs did not show over-representation within RISCs (S8B Fig). This lack of aggregate MHV68 miRNA enrichment in RISC, however, may be related to their non-canonical processing mechanisms [67]. Nevertheless, at the individual miRNA level, we found that the two most overrepresented miRNAs within RISC in the MHV68-infected cell line, HE2.1, were the MHV68 M1-2-3p and M1-6-5p miRNAs (S8B Fig, right), suggesting that at least some MHV68 miRNAs are overrepresented in RISC. Together, these findings show preferential loading of viral miRNAs into RISC across these three γHVs, representing a possible common strategy to exert strong influences on infected cell transcriptomes.

EBV miRNAs have high targeting efficacies
The disproportionately high association of viral miRNAs with RISCs could result from intrinsic properties of viral miRNAs that enhance loading into RISC. Nevertheless, high expression of viral miRNA and/or their targets could also contribute to the detection of greater numbers of miRNA-target RISCs. To account for the impact of miRNA and mRNA abundance on the level of AGO-bound miRNA-mRNA hybrids [68], we used an affinity constant calculation (dividing normalized hybrid counts by the normalized number of miRNAs and mRNA target transcripts Fig 3C)) as a quantitative metric for miRNA "targeting efficacy". We first assessed how accurately this targeting efficacy metric reflects intrinsic properties of miRNAs rather than cell specific factors. Unlike cellular miRNAs, nearly all viral 3'-UTR miRNA targets identified in Akata cells were also detected in SNU719 cells (Fig 3A), with the higher overall viral miRNA expression in SNU719 accounting for the additional miRNA targets identified in this system. Despite the overlap of specific targets, the relative abundance of each viral miRNA-mRNA hybrid detected in both cell lines was poorly correlated (Fig 3B), likely due to the differences in mRNA and/or miRNA expression in SNU719 and Akata cells. In contrast, viral miRNA targeting efficacies were strongly correlated between cell lines (Fig 3D). This finding also extended to cell miRNAs using miRNA/target interactions found in both cell lines ( Fig  3D). These data support our contention that this targeting efficacy metric is a quantitative measure of intrinsic properties of miRNAs and their target sites that influence RISC formation. We next applied targeting efficacy measures to assess the innate RISC loading properties of viral and cellular miRNAs in Akata and SNU719 cells. Across each biological replicate, the targeting efficacy distributions were substantially higher for EBV miRNAs than cell miRNAs in both SNU719 and Akata cells (Fig 3E and S2 Table). Extending this analysis to KSHV and MHV68 miRNAs, we similarly found that viral miRNAs have higher targeting efficacies than their cellular counterparts (S9A and S9B Fig and S2 Table). These results indicate that the higher proportion of viral miRNAs in RISCs relative to their expression levels is due to intrinsic properties of viral miRNAs that influence target selection and RISC loading. This suggests that viral miRNAs have evolved to be more effective inhibitors of their mRNA targets than host miRNAs.

EBV miRNAs target more accessible regions of mRNAs than cell miRNAs and form more thermodynamically stable hybrids
To assess biophysical properties underlying efficient viral miRNA RISC formation, we explored the nucleotide composition of miRNA/mRNA hybrids. Hybrid formation is guided in large part by the degree of complementarity between bases 2-8 of miRNAs and their targets, with loading into RISC being further boosted by the presence of an "A" opposite the first base of the miRNA [69]. These bases, collectively referred to as the "seed" region, display distinct complementarity patterns that are categorized into 9 different classes. With seed class being a known determinant of targeting effectiveness, we first tested whether there is a seed class bias for EBV miRNAs relative to cell miRNAs. As shown in S10 Fig, there is no discernable enrichment of EBV miRNAs in seed classes with higher known targeting capabilities. Further, we found that within each seed class, viral miRNAs have higher targeting efficacies (Fig 4A). These data indicate that differential seed class utilization does not explain the increased targeting efficacy of EBV miRNAs.
In addition to the extent of seed matching, miRNA binding is improved when target sites are flanked by A/Us [70]. Considering the 4 bases directly adjacent to the seed binding site (2 bases on each side of the seed binding region), we found that for both cell and viral miRNAs, sites with high A/U content (3+ flanking A/Us) were associated with improved targeting efficacy (Figs 4B and S11A). Nevertheless, EBV, KSHV and MHV68 miRNAs tended to bind sites with fewer flanking A/Us (Figs 4C and S11B), indicating that the number of flanking A/Us fails to explain the higher targeting efficacies of viral miRNAs.
Because secondary structure at target sites reduces accessibility and interferes with RISC formation [71], transcript regions with less intramolecular binding tendencies are typically more effective miRNA targets. We therefore assessed target site accessibility for each CLASH hybrid using the RNAplfold [72] subpackage of the Vienna RNA suite using a window size and maximum base pairing separation of 80 and 40 bases, respectively. This analysis revealed that EBV miRNA target sites are, on average, more accessible than the target sites of host miRNAs (Fig 4D). This indicates that EBV miRNAs evolved to target more accessible regions, providing one likely explanation for the observed higher targeting efficacies of EBV miRNAs.
The last feature that we assessed to determine the molecular basis of enhanced EBV miRNA targeting efficacies was predicted minimum free energies of miRNA-mRNA hybrid pairs (across the entire miRNA and its target). This analysis revealed that EBV miRNA-mRNA hybrids tend to form more thermodynamically stable interactions than their cellular counterparts (Fig 4E). Extending this analysis to KSHV and MHV68 miRNAs, we found that like EBV miRNAs, KSHV and MHV68 miRNAs generally form more thermodynamically stable interactions than their cellular counterparts (S11C Fig). Together, these analyses indicate that γHV miRNAs evolved with nucleotide sequence compositions that favor stronger hybrid interactions, providing another likely explanation for the enhanced targeting efficacy of viral miRNAs.

EBV miRNAs have a greater impact on their targets than their cellular miRNA counterparts
To test whether the higher overall targeting efficacies of EBV miRNAs translates to increased effectiveness in destabilizing their target mRNAs, we performed correlation analyses between the expression of each viral and cellular miRNA and their respective mRNA targets across EBV positive BL and GC datasets. Correlations were displayed in cumulative distribution plots (Fig 5). Shifts in distributions from permuted correlations of the same set of miRNAs with random mRNA targets were analyzed using the Kolmogorov-Smirnov (KS) test (Fig 5). In both BLs and GCs, viral miRNAs showed strong inverse correlations with their targets (BL, p < 2.1x10 -6 ; GC, p < 1.2x10 -5 ). In contrast, correlations between human miRNAs and their targets were less pronounced and not statistically significant. These results provide in vivo evidence of functional impacts of EBV miRNAs on targets identified by CLASH and they show that the higher targeting efficacies of EBV miRNAs likely translates into stronger functional influences on their targets.
While EBV miRNA targeting is generally spread across more transcripts, they tend to be directed towards transcripts regulated by fewer miRNAs (Fig 6C). This indicates that cellular miRNAs more often work in a cooperative fashion. The higher targeting efficacies of viral miRNAs may reduce their requirement for cooperative targeting, as they tend to target a repertoire of transcripts that are distinct from cellular miRNAs. As an extension of their enhanced targeting efficacy, viral miRNAs are not constrained by the necessity of cooperative targeting, allowing them each to serve a unique purpose, and increasing their influence over host cell transcriptomes.

EBV miRNAs interfere with immune signaling in the tumor microenvironment
Next, we sought to investigate the potential functional impact of the most influential EBV miRNA-target interactions in modulating the host cell environment. Using the targets of the top 20 EBV miRNA-target interactions as network seed genes, we applied enrichR[54] to query pathway enrichment for each seed gene and its 20 strongest interacting partners. As a control, we performed a similar analysis on the targets of the top 20 cell miRNA-target  [108]. Cumulative distribution plots were generated using ranked Spearman correlation coefficients for each species. As a control, expression levels of the miRNA of each hybrid pair analyzed was correlated with expression levels of a random mRNA across the same tumors. For each randomized interaction, 100 permutations were run. Correlation values of EBV miRNA-mRNA pairs were compared to randomized controls; P-values were generated using the KS test. https://doi.org/10.1371/journal.ppat.1009217.g005

PLOS PATHOGENS
interactions. This analysis revealed selective enrichment for Influenza and HIV infection, antigen processing and presentation (MHC class I), Ubiquitin & Proteasome degradation (which can impact antigen processing and presentation) and IFN-stimulated genes and ISG15 antiviral mechanisms as top EBV miRNA pathway hits (Fig 6D and S4 Table). This suggests that the EBV miRNA-target interactions within the top 20 highest targeting efficacy interface with adaptive and innate immune response pathways.
Previous studies have shown that EBV positive GCs have higher immune cell infiltration than their EBV negative counterparts [74], likely due to sporadic expression of lytic viral proteins within the tumor. Consistent with these studies, we used CIBERSORTx [75] to infer the immune cell infiltration in each tumor through deconvolution of immune cell signatures within the GC dataset (using only the microsatellite stable (MSS)-only cohort) and found higher levels of T-cell and macrophage subtypes in the EBV positive cohorts (S12A Fig). We also determined the diversity of infiltrating T-cells within these tumors by assessing the number of unique T-cell receptor sequences for each tumor. Utilizing MIXCR [76], we realigned each tumor RNA-seq dataset to all potential combinations of rearranged T-cell receptors and used VDJTOOLS [77] for QC and T-cell clone quantifications. This analysis showed that consistent with the finding of higher levels of T-cell infiltration in EBV positive GCs through CIBERSORTx analysis, we found a greater diversity of T-cell clones in EBV positive tumors than in EBV negative tumors (S12B Fig). Together, these results confirm the findings of The Cancer Genome Atlas Research Network [57], showing that EBV likely induces some level of adaptive immune response in EBV positive GCs.
While the presence of EBV clearly induces immune cell infiltration into EBV associated tumors, we hypothesized that the interactions of viral miRNAs with immune regulatory pathways in infected cells acts as a counter measure to mitigate immune cell-mediated targeting of virus-infected tumor cells. To functionally assess this possibility we first identified potential direct and indirect targets of viral miRNAs in the EBV+ BL and GC cohorts by correlating the sum of viral miRNA expression values with the expression levels of each cell gene (S5 Table). Gene Set Enrichment Analysis (GSEA) [78,79], using the pre-ranked BL and GC miRNA-gene correlations, revealed that viral miRNA expression negatively correlates with the IFNγ and TNFα signaling pathways (Fig 7A). The inverse relationship between IFNγ and TNFα signaling pathways and EBV miRNA expression suggests that EBV miRNAs mitigate IFNγ-and TNFα-mediated amplification of the immune response. To more directly assess the relationship between EBV miRNAs and adaptive immune response, we used CIBERSORTx [75] to infer immune cell infiltration in each tumor sample from the BL and GC datasets (S6 Table). This analysis showed that CD8 T-cells, CD4 T-cells, and macrophages inversely correlate with viral miRNA expression (Fig 7B). Assessing the infiltrating T-cell diversity across these tumors showed a strong inverse correlation between unique T-cell clonotype counts and viral miRNA expression (Fig 7C and S7 Table), further supporting a role of EBV miRNAs in interfering with the immune response to EBV infection within these tumors. These findings are most tightly associated with EBV miRNAs rather than other expressed EBV transcripts since EBV miRNAs show a stronger inverse relationship with T-cell clonotype numbers than viral pol II transcripts (S13 Fig). Together these analyses provide evidence that the targeting of components of immune response pathways by EBV miRNAs results in the dampening of innate and adaptive immune responses to viral infection in EBV positive cancers.

Discussion
By integrating miRNA and mRNA expression levels with hybrid quantifications obtained via CLASH, we used an affinity constant calculation to assess targeting efficacy for each interaction. Comparing targeting efficacies of viral and host miRNA-target interactions, we found that viral miRNAs generally bind their targets more effectively than host miRNAs. The elevated efficacy of EBV miRNA-mRNA interactions transcended seed match and the level of flanking AU content. Instead, we found that the complexes formed between EBV miRNAs and their targets were more thermodynamically stable and tended to occur at target sites with less secondary structure. In metazoans, broadly conserved miRNA target sites tend to facilitate the most effective targeting interactions. However, these sites exhibit a negative selection over time that tempers miRNA binding efficacy. Unlike targets of broadly expressed miRNAs, coexpressed targets of miRNAs that exhibit tissue specific expression tend to evolve effectivityreducing alterations in target sites [80], suggesting that the major role of miRNAs is to finetune gene expression rather than to significantly alter it [81]. The continued coevolution of cell miRNAs and their target genes presents an internal arms race that ultimately reduces the impact of the majority of miRNA-mRNA interactions. These constraints do not apply to viral miRNAs which have evolved to effectively target a large number of host transcripts, with some miRNAs such as BART7-3p effectively binding many target mRNAs, and others, such as BART3-3p, homing in on a single transcript (notably, this latter class of viral miRNA might

PLOS PATHOGENS
implicate singular targets that are critical to viral persistence and/or expansion in host cells, potentially presenting a vulnerability to the virus). The tendency of viral miRNAs to exhibit high targeting efficacies therefore arm them with the unmatched ability to influence the landscape of host RNA expression.
A longstanding challenge in deciphering miRNA function is the accurate assessment of its targets [82]. Because the most favorable interactions require only 7 bases of complementarity, even the most infrequent 3' UTR binding motifs are widely dispersed throughout the transcriptome. Still, the number of true miRNA binding sites is at least an order of magnitude lower than the number of seed complementary sequences in transcriptome 3' UTRs. The ability to discern true miRNA targets has improved dramatically with the advent and continued development of binding site prediction algorithms [70,82], and even more so with high throughput experimental approaches such as PAR-CLIP [83]. These tools have seeded studies of EBV miRNA functions, leading to the discovery of consequential mRNA targets [84] including mediators of apoptotic signaling, BBC3 (BART5-5p) [85] and BAD (BART20-5p) [86] and regulators of innate immunity, NDRG1 (multiple) [87], MICB (BART2-5p) [88], IFI30 (BART1 and BART2), LGMN (BART1 and BART2), and CTSB (BART1 and BART2) [89]. However, both putative and empirically identified viral miRNA targets have often been difficult to validate [90]. In our system, the BBC3 interaction was the 220 th most prevalent hybrid of BART5-5p (38 h.c.p.m.), CTSB was the 375 th most abundant BART1-5p target (6 h.c.p.m.), and combined, NDRG1 hybrids ranked 3681 st among all viral miRNA targets (25 h.c.p.m.). We were unable to detect any viral miRNA hybrids with IFI30, LGMN, or BAD, while MICB and IL12A were not sufficiently expressed in SNU719 cells. Overall, little overlap exists between targets identified among previously published high throughput EBV miRNA targetome studies [90], yet one consistently identified interaction is that between BART3-3p and IPO7, an interaction that has subsequently been used as a control in studies of EBV miRNA targeting [35,91]. In our study, we detect substantial levels of this interaction (6846 h.c.p.m.; ranked 1 st among BART3-3p targets), suggesting that the quantitative nature of this approach helps distinguish the targets of both viral and host miRNAs that are likely to be impacted the most.
A more general possible functional consequence of inhibiting ubiquitin ligase expression is interference with MHC class I antigen presentation. Previous studies have shown that EBV miRNAs collectively reduce host cell antigen presentation, culminating in the subversion of CD4+ and CD8+ T-cell surveillance [89,95,96]. Numerous viruses have evolved various mechanisms to interfere with this pathway. For example, the HIV NEF protein redirects HLA-A and HLA-B transportation [97], the Hepatitis B Virus X protein interferes with proteasomal activity [98], K3 and K5 of KSHV promote MHC-I endocytosis and destruction [99], MHV68 MK3 is a ubiquitin ligase that targets MHC-I for destruction [100], and both the U6 protein of HCMV and EBV encoded BNLF2A inhibit TAP-mediated peptide transport and subsequent MHC loading [101][102][103]. While viral protein expression is restricted in EBV-associated tumors, the presence of the virus still elicits an immune response and mutations that help drive tumors produce neo-epitopes that immune cells may recognize as foreign. The survival of EBV positive cells, including tumors, therefore, requires immune subversion or escape [104]. Previous studies have shown that EBV miRNA expression is elevated by inoculation of tissue culture grown EBV positive tumor cells in mice [105]. Furthermore, we found that EBV miRNAs inversely correlate with immune cell infiltrates in BLs and GCs, suggesting that an important role of the viral miRNAs is to interfere with immune surveillance, thereby helping facilitate the survival and expansion of the underlying tumor.
Altogether, our findings show that EBV miRNAs are uniquely effective inhibitors of target RNAs through facilitating high in vivo expression and the targeting of accessible sites with more favorable miRNA/target binding energies, thereby achieving conserved impacts on immune responses to infected tumor cell populations. Raw sequencing reads of latent and Zta-induced Akata cells (SRA accession: SRP042043) [109] were aligned to the combined hg38 + EBV reference transcriptome using kallisto [46]. The viral percentage of expressed mRNAs was calculated using the following equation,  [41,42]. CLASH sequencing reads were processed using the hyb bioinformatic pipeline [111] 2), we performed an α-AGO immunoprecipitation without crosslinking or ligation, followed by qRT-PCR for several selected microRNAs. The loading efficiency for each microRNA was calculated by taking the @Ct between AGO-IP and input sample (2 À ðaAgo Ct À Input Ct Þ ; x-axis). These values were compared to those obtained using CLASH and small fraction RNA-seq ( microRNA  Fig 4).