Next-Generation Sequencing Analysis Reveals Differential Expression Profiles of MiRNA-mRNA Target Pairs in KSHV-Infected Cells

Kaposi’s sarcoma associated herpesvirus (KSHV) causes several tumors, including primary effusion lymphoma (PEL) and Kaposi’s sarcoma (KS). Cellular and viral microRNAs (miRNAs) have been shown to play important roles in regulating gene expression. A better knowledge of the miRNA-mediated pathways affected by KSHV infection is therefore important for understanding viral infection and tumor pathogenesis. In this study, we used deep sequencing to analyze miRNA and cellular mRNA expression in a cell line with latent KSHV infection (SLKK) as compared to the uninfected SLK line. This approach revealed 153 differentially expressed human miRNAs, eight of which were independently confirmed by qRT-PCR. KSHV infection led to the dysregulation of ~15% of the human miRNA pool and most of these cellular miRNAs were down-regulated, including nearly all members of the 14q32 miRNA cluster, a genomic locus linked to cancer and that is deleted in a number of PEL cell lines. Furthermore, we identified 48 miRNAs that were associated with a total of 1,117 predicted or experimentally validated target mRNAs; of these mRNAs, a majority (73%) were inversely correlated to expression changes of their respective miRNAs, suggesting miRNA-mediated silencing mechanisms were involved in a number of these alterations. Several dysregulated miRNA-mRNA pairs may facilitate KSHV infection or tumor formation, such as up-regulated miR-708-5p, associated with a decrease in pro-apoptotic caspase-2 and leukemia inhibitory factor LIF, or down-regulated miR-409-5p, associated with an increase in the p53-inhibitor MDM2. Transfection of miRNA mimics provided further evidence that changes in miRNAs are driving some observed mRNA changes. Using filtered datasets, we also identified several canonical pathways that were significantly enriched in differentially expressed miRNA-mRNA pairs, such as the epithelial-to-mesenchymal transition and the interleukin-8 signaling pathways. Overall, our data provide a more detailed understanding of KSHV latency and guide further studies of the biological significance of these changes.


Introduction
Kaposi's sarcoma-associated herpesvirus (KSHV; also known as human herpesvirus 8) is an oncogenic human γ-herpesvirus involved in the pathogenesis of several AIDS-associated malignancies, including Kaposi's sarcoma (KS), primary effusion lymphoma (PEL) and multicentric Castleman's disease (MCD) [1][2][3][4]. As with other chronic viruses, infection with KSHV results in a number of changes in the host cells. Many of these are mediated by KSHV to evade the innate and adaptive immune responses, prevent cell cycle arrest, inhibit apoptosis, modulate cellular signaling pathways, and facilitate lifelong infection in the host. Conversely, other changes occur as a result of the host defense response to viral infection.
MicroRNAs (miRNAs) are a subtype of small (~22 nucleotides) noncoding RNAs found in all metazoan eukaryotes [5]. Genes encoding miRNAs are primarily transcribed by RNA polymerase II, generating imperfect stem-loop hairpin structures that are then processed by the cellular proteins Drosha, DGCR8, and Dicer to produce a miRNA duplex. This duplex is then incorporated into an RNA-induced silencing complex (RISC) that binds one miRNA strand to messenger RNA (mRNA) targets, usually in the 3' untranslated region (3' UTR) and induces translational inhibition or degradation (for a review, see [6]). While the changes in gene expression due to a given miRNA are often modest, they are nevertheless important post-transcriptional regulators because a single miRNA has the potential to alter an entire biological pathway by inhibiting numerous mRNAs at once through partial base complementarity [6].
The human genome encodes thousands of miRNAs [7]. A limited number of viruses, including KSHV and other herpesviruses, also encode their own miRNAs [8]. KSHV expresses 12 viral precursor miRNAs located within the latency-associated region and these yield 25 mature miRNAs [9,10]. Being non-immunogenic, virally encoded miRNAs are potentially attractive tools for viruses as they can interfere with the host machinery without being detected by the host immune system.
In addition, a number of KSHV-encoded proteins can effect substantial changes in host cell gene expression, either by directly acting on one or more steps in protein expression, or by indirect mechanisms (e.g. cell signaling pathways). In particular, KSHV proteins expressed during viral latency, such as latency-associated nuclear antigen (LANA) or viral FLICE inhibitory protein (vFLICE), can directly induce changes in expression of certain mRNAs or miRNAs to facilitate the latent infection. The changes in cellular miRNA expression can in turn affect the expression of target genes. In addition, changes in cellular miRNAs may occur as part of the host response to viral infection.
Transcriptome profiling techniques such as RNA-sequencing and microarrays have been used to study PEL cells, which are a valuable model system for KSHV infection [11]. PEL lines are largely latent but can have some lytic KSHV gene expression. However, a limitation of this approach is the lack of otherwise identical uninfected control lines. To bypass this restriction, differential expression studies using microarrays or high-throughput sequencing techniques have explored other KSHV infection models, such as KS lesion samples and normal skin tissues [12][13][14], mock and KSHV-infected primary lymphatic endothelial cells [15,16], or the KSderived epithelial SLK/SLKK (SLK+rKSHV.219) model [17,18]. Although previous studies have provided useful insights into KSHV infection and disease pathogenesis, to date there has been no report of a next-generation sequencing study combining both mRNA and miRNA expression profiles in KSHV-infected and uninfected control cells. Here, we present the outcome of a comprehensive analysis comprised of both small RNA-sequencing and poly-A enriched mRNA-sequencing in latently infected SLKK cells and uninfected control SLK cells. While these cells are derived from renal cancer and differ somewhat from the principal target cells in KSHV malignancies, advantages of this cell system are that SLKK cells can be compared to an uninfected counterpart and also that they have very low levels of spontaneous lytic KSHV replication. Thus we could focus on changes induced by latent infection, rather than those induced by an acute host response or expression of lytic KSHV genes. Furthermore, the combined analysis of miRNA and mRNA differential expression of SLK and SLKK cells enabled the identification of a number of viral-induced miRNA changes that were associated with expected inverse changes in the predicted target mRNAs, suggesting that miRNAs were mediating many of these changes. We further investigated the functions and pathways in which those miRNA-mRNA pairs were involved and highlighted their potential role in pathogen-influenced mechanisms. Overall, these findings provide evidence for a broad posttranscriptional regulatory network of miRNAs and their mRNA targets that will deepen our understanding of KSHV effects on the host cellular machinery and may inform further studies of specific miRNA-mRNA interactions in KSHV-infected cells.

Cell culture
Human KS-derived SLK and SLKK cells (also known as SLK+rKSHV.219) [19][20][21][22] were a gift from Dr. Don Ganem (UCSF, CA). They were expanded on receipt, frozen in liquid nitrogen, and stored in a cryogenic tank until used in the experiments described hereafter. Cells were thawed and maintained in Dulbecco's Modified Eagle medium supplemented with 10% v/v fetal bovine serum (Sigma-Aldrich, St Louis, MO) and 1% Penicillin/streptomycin/glutamine solution (Gibco, Carlsbad, CA). Additionally, KSHV-positive SLKK cells were periodically grown under selection with 10 μg/mL puromycin to maintain the viral episome. Not counting the time during which cells were frozen, SLK cells had been kept in culture for less than 3 months before being used for RNA isolation (RNA-Seq), transfection assays, or de novo KSHV infections. SLK cells infected with recombinant rKSHV.219 (SLKK cells) were originally selected over time using puromycin in order to reach a high latent viral expression. Therefore, again not counting the time that they were frozen, the SLKK cells were kept in culture for a longer period of time (6 months to a year) before being used for RNA isolation (RNA-Seq) and transfection assays. All experiments were done in triplicate independent cell cultures maintained at 37°C in humidified 5% CO 2 .
No. 4464051) and the double-stranded miRNA mirVana mimics for miR-708-5p, miR-409-5p and miR-409-3p (assay ID MC11161, MC13028, and MC12446, respectively; Ambion, Life Technologies) were transfected at a final concentration of 10nM using DharmaFECT 1 (Dharmacon, GE Healthcare, Lafayette, CO) according to manufacturer's protocol. Briefly, Dharma-FECT 1 and Opti-MEM (1:50 v/v ratio) were incubated at room temperature for 10 min. Separately, Opti-MEM and 10 pmol of miRNA mimics or miR-scramble were incubated for 10 min at room temperature. 400 μL of these two mixtures (1:1 v/v ratio) were then incubated for 20 min at room temperature and subsequently added to the cells with 1.6mL of DMEM containing 10% FBS. Cells were harvested for RNA 48 hrs after transfection.

RNA isolation
Total RNA was extracted from cells using miRVana miRNA isolation kit according to manufacturer's instructions (Ambion, Life Technologies, Carlsbad, CA). RNA abundance and integrity were determined after isolation using a Nanodrop-ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA), respectively. Only samples of total RNA with an RNA integrity number (RIN) >9 were further used for RNA-Sequencing and small RNA-Sequencing. All samples were stored at -80°C.

Preparation of microRNA libraries for deep sequencing
Small RNA libraries were constructed as previously described [23] using the TruSeq RNA Sample Prep kit (Illumina inc., San Diego, CA). Briefly, after running 2 μg of total RNA in a 15% urea-TBE gel (Invitrogen, Life technologies, Carlsbad, CA) for 1 hr at 200V, the 20 to 30 nucleotide RNA fraction was excised and eluted in 0.3 M NaCl. Following separation of the elute from the gel debris using a Spin-X-column (Thermo Fisher Scientific), the small RNA samples were precipitated in 100% ethanol and 1mg/mL glycogen, incubated at -80°C for 30 min, centrifuged at 14,000rpm for 25 min, washed with 75% ethanol, air dried and resuspended in RNAse-free water. Illumina TruSeq libraries were then prepared according to the manufacturer's protocol and the final RNA library concentration was measured by a Qubit 2.0 fluorometer using Qubit dcDNA HS Assay kit (Life Technologies). We verified the size of the products contained in the libraries using a high sensitivity DNA chip and an Agilent 2100 Bioanalyzer (Agilent Technologies). Finally, a total of six miRNA libraries (three SLK and three SLKK samples) were sequenced using the Illumina HiSeq platform.

Preparation of polyadenylated mRNA libraries for deep sequencing
Total RNA of samples used for small RNA-Sequencing was treated with Turbo DNA-free DNase I and Dynabeads (both from Ambion, Life Technology) to deplete samples from the residual DNA and to isolate the polyadenylated mRNA transcriptome, respectively. PolyA+ RNA libraries were then prepared with the ScriptSeq v2 RNA-Seq kit (Epicentre, Madison, WI). The final concentration and size distribution of the RNA libraries were measured by using a Nanodrop-ND-1000 spectrophotometer (Thermo Fisher Scientific), and by running a DNA 100 chip on an Agilent 2100 Bioanalyzer (Agilent Technologies). Finally, a total of six polyadenylated mRNA libraries (three SLK and three SLKK samples) were sequenced using the Illumina HiSeq platform.

Bioinformatic analysis
For small RNA sequencing, reads were aligned against known miRNAs from miRbase (version 19.0) using the software package SHRIMP [24] (http://compbio.cs.toronto.edu/shrimp/README) with miRNA specific settings (-o 1-n 2-r 30%-h 50%-local-Q-qv-offset 32-sam). To process paired-end sequencing, reads were aligned separately covering the mature miRNA both on the forward and reverse read and the obtained number of matches were averaged. The match counts were normalized and tested for differential expression using the edgeR package [25] with default settings. For mRNA sequencing, adapter sequences were trimmed using the Fastx toolkit. Reads were aligned against two target genomes, human (hg19) and KSHV (Genbank accession number NC_009333.1) using TopHat [26] to generate spliced alignments. Transcripts were assembled using Cufflinks and Cuffdiff [27] in order to reveal differentially expressed genes. To be more conservative and eliminate infinite fold differences in cases where there was a nil denominator, fold changes were calculated adding 0.01 read to both the numerator and denominator. Significant mRNA fold change was determined by an adjusted P-value lower than 0.05 based on the Benjamini and Hochberg multiple testing correction. In order to visualize miRNA-Seq profiles, sequencing reads were converted to bedGraph format using BEDtools. The output files were then uploaded and displayed using the UCSC Genome Browser (http://www.genome.ucsc.edu). R (http://www.Rproject.org) suite software was used for statistical analyses, heatmaps and scatter plots.

DNA isolation and PCR amplification
Genomic DNA was isolated from SLK and SLKK cells using the PureLink Genomic DNA Mini kit (Invitrogen) following the manufacturer's instruction. DNA concentration and integrity were determined after isolation using a Nanodrop-ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). All samples were stored at -20°C. Eight individual regions of the 14q32 locus were probed using a set of custom-designed primers at 10μM (S1 Table) and AmpliTaq Gold 360 Master Mix (Applied Biosystems). DNA amplifications were performed in a final volume of 25μL containing 100ng of DNA template. The cycling conditions were as follows: 10min at 95°C, 35 cycles of amplification at 95°C for 30 sec, 60°C for 30 sec and 72°C for 1 min, followed by 7 min at 72°C. The final products were electrophoresed on a 1% agarose gel, stained with ethidium bromide and visualize under UV light.
Pathway analysis and inverse correlation study of miRNA and mRNA expression levels Pathway analysis and inverse correlations between expression levels of differentially expressed miRNAs and their respective target mRNAs were analyzed using Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Redwood City, CA; www.ingenuity.com). This in silico analysis software reveals enrichment for molecular networks and signaling pathways. Moreover, we generated medium to high confidence miRNA target predictions and as well as experimentally observed miRNA-mRNA interactions using the IPA tool called "MicroRNA Target Filter" that integrates multiple target prediction algorithms such as TargetScan, TarBase, miRecords and the Ingenuity Knowledge Base. This allowed for the integration of miRNAs with mRNA targets, predicted and validated, that were both differentially expressed in SLKK vs. SLK cells. Opposite expression pairing between miRNA and mRNA levels was implemented to further refine the analysis. Further filtering options were applied such as the confidence parameter.
Gene set enrichment analysis (GSEA). GSEA [29] was performed by the JAVA program using the c3 collection of Molecular Signature Database, MSigDB (http://www.broadinstitute. org/gsea/msigdb/index.jsp). 221 miRNA families are each linked to a separate gene set consisting of putative targets containing a complementary 3'-UTR miRNA binding motif. The input gene list was generated by our mRNA-Seq differential expression analysis and was tested for enrichment in targets belonging to individual miRNA target genesets. A P-value and a false discovery rate (FDR), estimating the statistical significance of enrichment, were calculated for each miRNA gene set.

Distribution of cellular and viral miRNAs in KS-derived human SLKK cells by miRNA-Seq
Although KSHV only encodes for 25 mature miRNAs, they have been found to make up a significant proportion of the total mature miRNAs in PEL cells [11,[30][31][32]. We analyzed cellular and KSHV viral miRNA expression in latently infected SLKK cells by deep sequencing. We observed that a sizable fraction (~9%) of miRNA reads in SLKK cells were of KSHV origin (S2 and S3 Tables). The distribution of the different KSHV miRNAs in SLKK cells (Fig 1) was generally similar to that of PEL cells BC-1, BC-3 and BCBL-1 reported previously [11]. For example, in SLKK cells, KSHV-miR-K12-4-3p and 8-3p comprised a similar proportion of KSHV miRNAs as in PEL cells (~26% and~13%, respectively). However, there were some differences; KSHV-miR-K12-6-3p (2%) was expressed at a lower level, and KSHV-miR-K12-10a-3p (40.5%) was expressed at a higher level in SLKK cells compared with PEL cells (20-43% and <1%, respectively). It should be noted that substantial differences in KSHV miRNA distribution were previously observed even within the same cell type (B-cells, i.e. BC-1, BC-3 and BCBL-1) [11,30]. It is therefore not unexpected to observe similarities but also dissimilarities between the SLKK viral miRNA profile and that of PEL cells. Such differences may also arise from cell type differences between B-cells and SLK cells, the use of rKSHV.219 in the SLKK cells, or linker ligation biases.

Comparison of expression profiles reveals more down-regulated than up-regulated cellular miRNAs in KSHV-infected cells
We compared the cellular miRNA expression levels between the small RNA libraries of noninfected and KSHV-infected cell cultures. For the purpose of this analysis, we applied thresholds of 2 times fold change (FC) in expression (-1 and 1 log 2 FC) and a P-value <0.05. Overall, we found 153 human mature miRNAs significantly altered by KSHV infection in SLKK vs. SLK cells (Fig 2A and S1 Fig). Interestingly, we observed substantially more suppressed than induced miRNAs, 111 vs. 42 respectively. In terms of expression fold change, this tendency was also apparent amongst the miRNAs with the greatest changes, which extended to an 11-fold down-regulation but only to a 4-fold up-regulation in log 2 scale (Fig 2A). Amongst the most dysregulated miRNAs that were expressed in relatively high amounts (threshold of mean miR count > 10), miR-224-5p (-11 log 2 FC), miR-452-5p and miR-887 (both -9 log 2 FC), were significantly down-regulated, whereas the greatest relative induction was seen for miR-708-5p/3p and 3614-5p (all~2 log 2 FC) ( Fig 2C).

The 14q32 cluster is down-regulated by KSHV infection
It is noteworthy that 42 of the 111 miRNAs that were observed to be repressed by KSHV infection are located in the same genomic region, the 14q32 miRNA cluster (S2 Fig). For example, 14q32 miRNAs, miR-410, miR-409-5p and -3p, were significantly repressed -6 log 2 FC upon KSHV infection (Fig 2C). Of note, while both arms of the miR-409 stem-loop were downregulated, miR-409-3p had a higher relative abundance than miR-409-5p. Interestingly, there has been a previous report of frequent loss of 14q32 in PEL cell lines [33], which suggests that down-regulation of this region at the DNA level may be beneficial for the survival of KSHVinfected cells. To test whether the down-regulation in SLKK cells was similarly due to the loss of DNA encoding these miRNAs, we interrogated eight individual regions of the genomic DNA of SLK and SLKK cells by performing PCR (S3 Fig). We found that all tested portions of the 14q32 cluster were present in uninfected as well as KSHV-infected cells, suggesting posttranscriptional mechanisms lead to the substantial down-regulation of 14q32 miRNAs.

mRNA-Seq analysis of SLKK cells
In parallel with these deep sequencing based miRNA studies, we also investigated the global changes in gene expression associated with chronic latent KSHV infection in order to assess the effects of the miRNA changes on their predicted targets. We compared the cellular mRNA expression levels of the poly-A enriched RNA libraries between non-infected SLK and infected SLKK cells. In all, 1,570 cellular genes,~5% of all expressed genes, were significantly altered with either -1 or 1 log 2 FC (2 times linear FC) and an adjusted P-value (FDR) <0.05 ( Fig  2B and S4 Fig). This included 980 repressed and 590 up-regulated mRNAs (Fig 2B). Some of these differentially expressed genes have been previously reported to be associated with KSHV infection. For example, in concordance with previous findings, DUSP1, JUN, and PLAUR were repressed whereas MDM2, CD55, CXCL1, IL6ST and ITGB4 were induced by KSHV infection [34][35][36][37]. Also, KSHV infection induced VCAM-1 (3 log 2 FC), a vascular endothelial marker that has been linked to tumor angiogenesis and oncogenesis [38].
Some of the mRNAs that were decreased upon KSHV infection in SLKK cells have been previously described as targets for KSHV-encoded viral miRNAs. For example, the most abundant KSHV miRNA in SLKK cells, miR-K12-10a, is known to target TWEAKR (TNFRSF12A) [34], which was significantly repressed in SLKK vs. SLK cells (-2.0 log 2 FC and FDR = 0.001; Fig 2B). Another KSHV miRNA, miR-K12-11, down-regulates a translational regulator of interleukin 6, C/BEPβ [39], which was also decreased in SLKK compared to SLK cells (-1.3 log 2 FC and FDR = 0.016). It is quite plausible that the decreases in these target genes are the result of KSHV miRNA expression.

Validation of miRNA differential expression patterns by Taqman assay
To validate the results observed in our miRNA sequencing data, we carried out Taqman assays (qRT-PCR) on several human miRNAs, with a focus on those with the greatest fold changes. Amongst the most significantly down-regulated miRNAs with high read counts were miR-224-5p, miR-410, miR-409-3p and -5p, whilst the most up-regulated miRNAs were miR-3614-5p, miR-708-3p and -5p. Because of the prominence of the miRNA repression, we also wanted to verify the differential expression of a miRNA that was less strikingly down-regulated, such as miR-212-3p. We also examined the change in miR-210 by qRT-PCR. This miRNA is of particular interest because it is induced by hypoxia-inducible factor (HIF) and because KSHV infection is known to increase HIF levels [40]. The levels of miR-210 by miRNA-Sequencing were slightly increased by KSHV infection, although the change did not achieve significance (Pvalue >0.05). Fig 3A shows the log 2 -transformed expression fold change of these miRNAs between SLKK and SLK cells, assessed by both miRNA-Sequencing and Taqman assay. The differential expression of all nine miRNAs tested by Taqman assay (eight of which met the RNAseq P-value cut-off) was significant and consistent with our sequencing data. Pearson correlation coefficient (r) between the two measurement methods of the eight selected miRNAs indicates strong and significant positive correlation (r = 0.98, p < 10 -5 ).
The relative expression of the 4 most down-regulated miRNAs seen in latently infected SLKK cells, 3 of which belong to the 14q32 cluster, was also examined following de novo KSHV infection of SLK cells. At 5 days post infection (5dpi), we found that 2 of these 4 miRNAs (miR-224-5p and miR-409-3p) were similarly down-regulated (n = 3) (Fig 3B). Also, miR-409-5p showed a trend for down-regulation, while miR-410 did not significantly change. A lower KSHV episome copy number in the de novo infected SLK cells as compared to SLKK cells could account for these changes. These results provide additional evidence that a proportion of the miRNA changes in SLKK vs. SLK cells are in response to viral infection. These results also suggest that some of the changes may only be observed in cells latently infected for a period of time.
Moreover, miR-708-5p, which was found up-regulated~3 log 2 FC by miRNA sequencing and~2 log 2 FC by qRT-PCR, is known to target pro-apoptotic protease caspase-2 (CASP2)  A. miRNA expression in latently infected SLKK cells compared to uninfected SLK cells and the mRNA expression of their corresponding targets (Taqman and qRT-PCR assays). miRNAs and mRNAs are shown in light and dark grey, respectively. miR-409-3p targets fibrinogen beta (FGB) and radixin (RDX). miR-409-5p is predicted to target MDM2. miR-708-5p targets caspase-2 (CASP2) and is predicted to target leukemia inhibitor factor (LIF). P-values were calculated using Student t-test. *, ** and *** indicate P <0.05, 0.01 and 0.001, respectively. B. Effect of miRNA transfection on target regulation. Scramble control (miR-scramble) and miRNA mimics for miR-409-3p, miR-409-5p and miR-708-5p were transfected in either SLK or SLKK cells. Quantitative real-time polymerase chain reaction revealed the target expression for these miRNAs and the relative value is presented as the mean ± standard deviation based on 3 independent experiments. Statistical analysis and its representation are as in Fig 4A. doi:10.1371/journal.pone.0126439.g004 [43]. Consistent with this, the qRT-PCR data demonstrated a small but significant reduction in CASP2 levels in KSHV-infected SLKK cells compared to control SLK cells. According to Tar-getScan, miR-708-5p is also predicted to bind the 3'UTR of leukemia inhibitory factor LIF mRNA and indeed we found it to be significantly down-regulated by qRT-PCR (Fig 4A). The strong reproducibility between both our miRNA-Seq and mRNA-Seq, and our qRT-PCR analyses provides evidence that the RNA-Seq was very robustly analyzed and enabled a valid assessment of the miRNA and gene expression profiles.
In order to further investigate whether the changes in mRNA targets in SLKK cells could be a result of differential miRNA expression, we carried out transfection assays to re-express miR-NAs that were decreased in either SLKK cells (miR-409-3p and -5p) or SLK cells (mir-708-5p), as compared to their counterparts. Transfection of miR-409-3p in SLKK cells led to a decrease in two of its known targets, fibrinogen beta and radixin. Additionally, transfection of miR-409-5p in SLKK cells led to a significant decrease in the expression of its predicted target, the p53-inhibitor MDM2. This is the first report providing experimental evidence that miR-409-5p targets MDM2. Finally, transfection of miR-708-5p in SLK cells resulted in a dramatic (80%) decrease in its target, caspase-2; it also led to the significant decrease of predicted target LIF (~40%). This provides evidence that the up-regulation of miR-708-5p in SLKK cells is likely causing the observed down-regulation of caspase-2 ( Fig 4B).

Two integrated analyses of miRNA and mRNA show post-transcriptional regulatory networks relevant to KSHV pathogenesis
Based on the above results, we hypothesized that at least some of the changes in mRNA expression levels in infected cells are caused by the differentially expressed miRNAs. We further hypothesized that these miRNA-mRNA pairs may play a role in KSHV infection or in the cellular response to chronic virus infection, and also in the pathogenesis of KSHV-related diseases. To explore this, we used MSigDB, a function of Gene Set Enrichment Analysis (GSEA), to investigate the enrichment of the differentially expressed genes in targets belonging to distinct miRNA families. A miRNA family is determined by its seed region, a 6-8 nucleotide sequence that is located at the 5' end of the miRNA and is able to bind target mRNA, promoting its degradation. Using our RNA-Seq data, we identified both up-regulated and down-regulated genes in the SLK/SLKK model. We then compared the GSEA output miRNA list to our miRNA-Seq data and selected for those miRNAs whose up-regulation is associated with the observed target down-regulation and then those whose down-regulation was associated with target up-regulation. This analysis allowed us to pair 4 up-regulated miRNAs with 74 down-regulated targets and 17 down-regulated miR-NAs with 195 up-regulated targets overall (Table 1). These findings support the hypothesis that some of the changes in gene expression from cells latently infected with KSHV are caused by reciprocal expression changes in the miRNA that target those genes.
A complementary approach for integration of miRNA and mRNA datasets is to use Ingenuity Pathway Analysis (IPA), which has an extensive collection of both known and predicted targets. IPA gathers comprehensive miRNA targeting information by combining multiple target prediction tools (TargetScan, miRecord and TarBase). Fig 5 summarizes the IPA data mining process using both sequencing datasets. Because miRNAs are known to effect relatively small changes in target mRNA expression and because changes induced by poorly expressed miR-NAs might be particularly hard to detect at the target level, we focused on the medium to highly expressed miRNAs that were significantly affected by KSHV infection.
Of the 153 differentially expressed miRNAs, 89 met the read count cut-off of 60 reads across all replicates (i.e. an average of 10 per replicate). IPA then identified 50 of these 89 miRNAs that had either predicted or validated targets (a total of 12,964 targets). Of these 50 miRNAs, 48 had targets in the list of differentially expressed (2 times linear FC and FDR <0.05) mRNAs generated by our sequencing analysis, for a total of 1,117 targets. Interestingly, each of the 48 miRNAs targeted at least one gene that was affected in the expected direction. In total, the differential expression of 73% of those target genes (814 out of 1,117) was inversely correlated to that of their respective miRNAs, suggesting that at least a portion of the differential mRNA expression was a result of changes in miRNA expression.
To adopt a conservative approach, we further filtered our dataset using experimentally observed or highly predicted target correlations. After this filter, we still had 45 miRNAs targeting 480 mRNAs, with each paired interaction being inversely correlated (Table 2 and Fig 5). Of these, 8 miRNAs were paired with 34 target genes that had been experimentally determined (Table 3), and 17 miRNAs had targets directly involved in pathogen-influenced signaling pathways (Fig 6). We calculated that a third of IPA miRNA filter results overlapped with those of MSigDB. Also, two third of the MSigDB matched IPA analysis, confirming the effectiveness of the IPA analysis.
Using these 480 differentially expressed mRNAs as input, pathway analysis revealed an enrichment in several signaling pathways that could be relevant to KSHV pathogenesis, such as "regulation of the epithelial-to-mesenchymal transition pathway", "IL-8 signaling", Using our list of differentially expressed mRNAs, we predicted which miRNAs would be upstream inhibitors with the Molecular Signatures Pathways Database (MSigDB). The overlap between the MSigDB output and our miRNA-Seq data highlighted 21 miRNAs that were differentially expressed in SLKK compared to SLK cells. Columns identify their names, the target sequence used to identify base complementarity with the input mRNA, the total number of putative target genes for a given miRNA, the number of putative targets that overlap with our input mRNA list, the percentage ratio between the overlap and the total number of putative targets, and its associated P-value and false discovery rate (FDR), correcting for multiple hypothesis testing. doi:10.1371/journal.pone.0126439.t001 "clathrin-mediated endocytosis signaling" and "Wnt/β-catenin signaling" (Table 4). Of note, EMT plays a central role in cancer progression and can be promoted by KSHV through Notchinduced ZEB1, a transcription factor up-regulated in our dataset [45]. Also, the Wnt/-catenin pathway, which regulates tumor development and cell proliferation, is altered by latencyassociated nuclear antigen LANA, a major KSHV latency protein suggesting multiple viral mechanisms that achieve the same result [46]. Finally, the enrichment of the interleukin-8 signaling pathway is noteworthy in that KSHV ORF74, which plays a key role in KS pathogenesis, is a constitutively active viral homologue to the interleukin-8 receptor. Altogether, these findings further support the role of miRNAs and their transcript targets in KSHV pathogenesis and/or the response to cells to chronic viral infection.  Using our lists of differentially expressed miRNAs and mRNAs as input for the Ingenuity Pathway Analysis (IPA), we found that 45 differentially expressed miRNAs target 480 mRNAs (medium and high confidence). Only mRNAs and miRNAs with opposite differential expression are shown in this table.
Columns identify the miRNA name, its log 2 -transformed expression fold change between SLKK and SLK cells, the number of identified targets, and the five or less most differentially expressed targets. In bold are miRNA prediction overlapping with GSEA analysis (Table 1). doi:10.1371/journal.pone.0126439.t002

Discussion
In this study, we used next-generation sequencing to survey both mRNA and miRNA expression profiles in SLKK cells latently infected with KSHV as compared with the uninfected SLK cells. This allowed us to determine the degree to which changes in gene expression profiles could be attributed to changes in miRNA expression profiles. Looking at the broad trends in differential miRNA profiles, the majority of miRNAs that met our selection criteria (2 linear FC and P0.05) were down-regulated in the KSHV-infected cells. This is consistent with previous reports that lymphatic endothelial cells infected with KSHV and primary B-cells infected with Epstein-Barr virus (EBV) show a substantially higher number of down-regulated than upregulated miRNAs [47,48]. One possible explanation was that KSHV-encoded miRNAs were produced at such a high level that they competitively interfered with the processing of the cellular miRNAs. However, in the SLKK cells, viral miRNAs made up only~9% of the total miRNA fraction, suggesting that this was not the principal mechanism. Considering alternative explanations, other viruses have been shown to globally reduce cellular miRNA levels by the saturation of exportin-5 (XPO5), which is involved in transit of pre-miRNA to the cytoplasm [49], and it is possible that this may be occurring in KSHV-infected cells.
An unexpected finding was that more than a third of the significantly down-regulated miR-NAs were located within 50kb of each other at the 14q32 miRNA cluster, which has previously been linked to cancer development. In fact, some of the 14q32 miRNAs that were easily detectable in the uninfected cells were undetectable in the KSHV-infected cells. We further validated the down-regulation of 14q32 miRNAs miR-409-3p, miR-409-5p, and miR-410 by qRT-PCR. Table 3. miRNA-mRNA pairs found inversely correlated and validated in the literature.

miRNA
Known Target  Reference   Up-regulated miRNAs/Down-regulated targets   let-7i-5p  ATAD3B, BCL2L1, CCND1, CDIPT, DUSP23, FAM105A, F2,  GAK, GPR56, GYS1, MARS2, POLD2, POM121/POM121C A substantial down-regulation of 14q32 miRNAs has similarly been observed in PEL cell lines including BCBL-1, BC-1 and BC-3 [11]. Because of the extent of the down-regulation in the current study, we wondered whether a deletion of this region at the DNA level might have occurred, as it has been reported in a number of PEL cell lines and primary cases of PEL [33]. In SLKK cells however, DNA was present in each of eight regions that we probed, suggesting that repression of 14q32 miRNAs is largely due to epigenetic silencing at the DNA level or possibly post-transcriptional mechanisms. Also, the observation that 14q32 miRNAs are frequently down-regulated in the setting of KSHV infection and that it can occur either through chromosomal deletions or through other mechanisms suggests that down-regulation of these miRNAs are advantageous for survival of KSHV-infected cells.
A strength of the study is that we have simultaneously used deep sequencing to survey paired changes in miRNA and mRNA targets, in both infected and infected cell lines without transfection or other manipulation. However, a weakness of this approach is that miRNA changes are not the only factors that can modulate mRNA in this system and even where paired changes are seen, this does not prove causality; KSHV-encoded proteins, as well as manifestations of the host response, can also affect cellular mRNA through mechanisms other than the miRNA silencing machinery. Nonetheless, our analysis provides evidence to suggest that a substantial proportion of the observed mRNA changes were caused by miRNA dysregulation. We compiled an extensive list of strongly correlated miRNA-mRNA pairs. Of 48 miRNAs whose level met our cut-off criteria for modulation, there were predicted or experimentally determined targets that were also modulated in the SLKK cells, for a total of 1,117 targets. Of these, 73% (814 targets) were modulated in the direction predicted by the miRNA change. These modulated mRNAs include 36 mRNAs that had been previously reported as known targets for 9 miRNAs dysregulated in SLKK vs. SLK cells. We further transfected 3 miRNAs that were decreased in SLK or SLKK cells and examined changes in 5 respective mRNA targets, including 2 novel mRNAs (LIF for miR-708-5p and MDM2 for miR-409-5p). In each case, transfecting the miRNA resulted in a predicted decrease in the target gene(s). Moreover, based on our current understanding, changes in some of these mRNAs would benefit KSHV infection. For example, a decrease in caspase-2, a proapoptotic miR-708-5p target, would benefit maintenance of KSHV latent infection [44]. Taken together, these results are thus consistent with the hypothesis that changes in miRNAs are driving at least some of these mRNA changes. Furthermore, these results provide some additional confidence in the use of paired miRNA-Seq and mRNA-Seq in the SLK/SLKK system to identify mRNA targets that may be relevant to KSHV pathogenesis. However, it is unlikely that all of the changes in any given mRNA target are directly and uniquely caused by a specific miRNA, as biological systems are often redundant.
It is worth keeping in mind that while miRNA dysregulation could result from KSHV infection itself, it may also be due to a host response to infection, or external factors such as genetic drift in the established infected cell lines. Interestingly there was relatively little overlap between the miRNA changes observed here in latently-infected cells and those reported to occur in the first few hours after de novo KSHV infection, which would be expected to predominantly reflect the antiviral innate immune response or KSHV genes expressed early in infection [16]. This suggests that our findings are largely related to the effects of long-term KSHV latency. Looking at the miRNAs overexpressed here in SLKK cells, previous reports provide evidence linking certain miRNA changes to a direct effect of latent KSHV infection. For example, The 480 differentially expressed genes targeted by miRNAs between SLKK and SLK cells were used as input for the Ingenuity pathway analysis (IPA).
Here, we highlighted 4 of the top 10 enriched pathways that are relevant to KSHV pathogenesis. Columns identify the pathway name, its associated Pvalue, the ratio between dysregulated genes and genes in the pathway, which dysregulated target gene and miRNA are involved. Genes or miRNAs in bold and normal font were found up-regulated and down-regulated, respectively, by chronic KSHV infection.
doi:10.1371/journal.pone.0126439.t004 miR-146b has been shown to be induced by STAT3 [50], which is phosphorylated in cells with latent KSHV infection [51]. miR-146b in turn can down-regulate Toll-like receptors and other aspects of the immune response to pathogens and it may thus help KSHV evade immunity [52]. Moreover, a STAT3 knock-down study showed that STAT3 induces miR-21-3p and miR-29c-5p, which were found up-regulated in SLKK cells; and STAT3 represses miR-2276, which was found down-regulated in our study [53]. Also, while the change in miR-210 was not sufficient to meet our original P-value miRNA-Seq threshold, it was nonetheless found significantly up-regulated based on the Taqman assay. KSHV LANA has been shown to up-regulate hypoxia-inducible factor (HIF) [54], a direct inducer of miR-210 [55,56], so the increase in miR-210 was predicted based on this mechanism. How various cellular miRNAs are controlled by KSHV and the potential role that KSHV proteins may have at the post-transcriptional level are important questions that remain to be explored in the future. Looking further at individual miRNAs dysregulated in SLKK cells and their target mRNAs, several stand out as potentially important in KSHV infection. For example, compared to SLK cells, SLKK cells manifested an increase in miR-708-5p along with a significant down-regulation of pro-apoptotic protease caspase-2 (CASP2). Through its induction of apoptosis, CASP2 can have an important anti-viral effect, and inhibition of this activity can thus promote KSHV infection. Also, miRNAs miR-485-5p and miR-409-5p are down-regulated in KSHV-infected cells and this directly correlates with the up-regulation of their predicted target, MDM2, which mediates the degradation of p53 [57,58], an inducer of apoptosis. By this mechanism, MDM2 may suppress virus-induced apoptosis, thereby increasing the survival of KSHV-infected cells. Two other down-regulated miRNAs, miR-299-3p and miR-381-3p, are predicted to suppress expression of integrin V (ITGAV) ( Table 2 and Fig 6), which was in turn up-regulated in SLKK cells. ITGAV has been shown to enhance the ability of KSHV-infected cells to infect new cells through binding to the extracellular matrix [59]. Interestingly, miR-409-5p, miR-299-3p, and miR-381-3p are all located in the 14q32 cluster that is broadly down-regulated in SLKK cells as well as PEL cells [33]. Additional studies will be needed to clarify the biological effects of these various miRNA changes.
In this study, we also used pathway analysis tools to identify some of the most affected pathways. Overall, we found 45 miRNAs and 480 paired mRNA targets that were either experimentally determined or highly predicted and whose changes inversely correlated with the mRNA changes. Of these, 17 miRNAs and 28 putative mRNA targets have been previously identified as being involved in pathogen-related signaling (Fig 6). We also noted several fundamental pathways that were significantly enriched in the 480 mRNA targets. These include the regulation of the epithelial-to-mesenchymal transition (EMT) pathway, IL-8 signaling, clathrin-mediated endocytosis signaling and Wnt/-catenin signaling (Table 4).
We used SLK and SLKK cells in the current studies because they provided the ability to compare KSHV-infected and uninfected cells and because of the tight control of KSHV latency in these cells. Experimentally derived from a KS biopsy, the uninfected SLK cell line was initially thought to have endothelial features matching those of KS and it was anticipated to become a KS cell model since it can be stably and latently infected by KSHV [21,22]. However, it was later shown to be indistinguishable from Caki-1, a clear-cell renal-cell carcinoma line of epithelial origin [20]. It is worth noting that SLK cells stably infected with recombinant rKSHV.219 (SLKK) were selected over time to only consist of KSHV-infected cells and some of the changes in SLKK cells as compared to SLK cells could be the result of genetic drift. However, some of the miRNA changes observed in the SLKK cell line were also observed upon de novo infection and therefore genetic drift could not account for all the changes. Thus, the SLK/SLKK system remains a valuable model to study host effects from latent KSHV infection and a number of changes observed here have also been observed in other studies of miRNA or mRNA expression in KSHV-infected cells [18]. This said, a number of complex factors may be contributing to specific observed changes, and additional studies will be needed to tease out these relationships.
In summary, we used deep sequencing to examine the small RNAs and mRNAs expressed by SLKK and SLK cells, and we identified a number of miRNA-target mRNA pairs that play a significant role in KSHV pathogenesis and in some cases, KSHV-induced diseases. We also observed that the majority of cellular miRNAs that are significantly altered by KSHV are downregulated, including most members of the 14q32 miRNA cluster, a genomic locus linked to cancer. To our knowledge, this is the first comparative study correlating changes at the miRNA and mRNA level between cells with latent KSHV infection and uninfected cells. This survey can help inform additional studies to better understand the mechanisms responsible for the specific changes observed and their role in the establishment of KSHV infection, the host response to this infection, and the pathogenesis of KSHV-induced tumors.
Supporting Information S1 Fig. Expression of miRNAs in the analyzed small RNA libraries. A. Scatterplot of the transcriptome of KSHV-infected versus uninfected cells. The X-axis represents the average of the normalized miR read counts across all samples. The Y-axis represents the logtransformed expression fold change between SLKK and SLK cells. Blue dots are differentially expressed genes with a P <0.05, while grey dots are those with a P 0.05. B. Heatmap of miRNA changes between SLK and SLKK cells in the six individual replicates (three SLK and three SLKK). Using uncentered Pearson correlation as the distance metric, we generated an unsupervised hierarchical heatmap where each row represents a miRNA and each column represents a biological replicate. Red, black and green denote low, median and high relative miRNA expression, respectively. The X-axis represents the average of the normalized read counts across all samples. The Y-axis represents the log 2 -transformed expression fold change between SLKK and SLK cells. Blue dots are differentially expressed genes with a FDR <0.05, while grey dots are those with a FDR 0.05. B. Heatmap of mRNA changes between SLK and SLKK cells in the six individual replicates (three SLK and three SLKK). Using uncentered Pearson correlation as the distance metric, we generated an unsupervised hierarchical heatmap where each row represents a mRNA and each column represents a biological replicate. Red, black and green denote low, median and high relative miRNA expression, respectively. (TIFF) S1 Table. Primer sequences for PCR amplification of the 14q32 genomic region. Eight individual regions (~450bp) of the genomic DNA at 14q32 were probed using these custom-designed primers. Columns identify the primer name, direction, sequence as well as amplicon size. (DOCX) S2 Table. KSHV miRNA distribution in SLKK cell line. This table shows the read count and percentage representation of 25 mature KSHV miRNAs expressed in SLKK cells. (DOCX) S3 Table. Reads statistics for miRNAs expressed in KSHV-positive SLKK cells. Three independent experiments are displayed (A, B and C). Columns identify the replicate number, raw read counts for human or mature miRNAs, number of mapped miRNAs per species, which percentage these KSHV reads represented for each replicate or overall. These counts were obtained by mapping reads to either KSHV mature miRNAs or human mature miRNAs. (DOCX)