A Non-Coding RNA Network Involved in KSHV Tumorigenesis

Regulatory pathways involving non-coding RNAs (ncRNAs), such as microRNAs (miRNAs) and long non-coding RNAs (lncRNA), have gained great relevance due to their role in the control of gene expression modulation. Using RNA sequencing of KSHV Bac36 transfected mouse endothelial cells (mECK36) and tumors, we have analyzed the host and viral transcriptome to uncover the role lncRNA-miRNA-mRNA driven networks in KSHV tumorigenesis. The integration of the differentially expressed ncRNAs, with an exhaustive computational analysis of their experimentally supported targets, led us to dissect complex networks integrated by the cancer-related lncRNAs Malat1, Neat1, H19, Meg3, and their associated miRNA-target pairs. These networks would modulate pathways related to KSHV pathogenesis, such as viral carcinogenesis, p53 signaling, RNA surveillance, and cell cycle control. Finally, the ncRNA-mRNA analysis allowed us to develop signatures that can be used to an appropriate identification of druggable gene or networks defining relevant AIDS-KS therapeutic targets.


INTRODUCTION
Non-coding RNAs (ncRNAs) are RNA transcripts that do not encode proteins and based on the length can be divided into two classes: small ncRNAs (sncRNAs), with transcripts shorter than 200 nucleotides, and long ncRNAs (lncRNAs), with transcripts longer than 200 nucleotides (1). Regulatory pathways involving ncRNAs, such as microRNAs (miRNAs), belonging to the class of sncRNA, and lncRNAs have gained great relevance due to their role in the control of gene (mRNA) expression. Different modes of interactions between lncRNAs and miRNAs have been reported: miRNA decay of lncRNAs, lncRNAs competing with mRNAs to bind to miRNAs, lncRNAs competing with miRNAs to bind to mRNA, and lncRNAs being shortened to miRNAs (2,3). All these interactions regulate the expression levels of mRNAs and in turn affect core protein signals, resulting in changes in the physiological functions of cells.
Kaposi's sarcoma (KS) is an AIDS-associated malignancy caused by the KS herpesvirus (KSHV). Despite the reduction of its incidence since the implementation of anti-retroviral therapy (ART), KS continues to be a global, difficult-to-treat health problem, in particular for ART-resistant forms (4,5). KS is characterized by the proliferation of KSHV-infected spindle cells and profuse angiogenesis (6).
The life cycle of KSHV has two well-defined phases: latent and lytic. In the latent phase, the virus expresses a few genes involved in viral persistence and host immune evasion. During the lytic phase, which is triggered by environmental and/or physiological stimuli, the viral genome replicates and new virions are formed (7). At this stage, KSHV is particularly effective at exploiting host gene expression for its own benefit. In this sense, the coevolution of the virus and its host has developed an intricate association between the virus genome, with its coding genes and non-coding genes, and the host RNA biosynthesis machinery (8). To the point that KSHV can seize control of RNA surveillance pathways, such as DNA damage response (DDR), pre-mRNA control machinery and the Nonsense-mediated mRNA decay (NMD), to fine-tuning the global gene expression environment throughout both phases of infection (7,9,10).
A recent study of KSHV-infected TIVE cells using wild-type and miRNA-deleted KSHV in conjunction with microarray technology to profile lncRNA expression found that KSHV can deregulate hundreds of host lncRNAs. These data established that KSHV de-regulates lncRNA in a miRNA-dependent fashion (11).
Using deep RNA sequencing of KSHV Bac36 transfected mouse endothelial cells (mECK36) and tumors (12), we have previously analyzed the host and viral transcriptome to characterize mechanisms of KSHV-dependent and -independent sarcomagenesis, as well as the contribution of host mutations (13). We now decided to study in this model, in a genome-wide fashion, the ncRNAs landscape to better understand the relationship between mRNAs, lncRNAs, and miRNAs in shaping KSHV tumorigenesis mechanisms.
This study allowed us to identify the most relevant host lncRNAs involved in KSHV tumorigenesis through the mouse KS-model (Malat1, Neat1, H19, and Meg3). In addition to having common target genes, pathway analysis showed that the four lncRNAs also share common related processes, mainly associated with cancer and viral infections, which would contribute with a network of genepathways closely associated with KSHV oncogenesis. We also showed evidence of the most frequent viral lncRNAs transcripts expressed in our model.
On the other hand, small RNA-sequencing and miRNA analysis revealed a high proportion of upregulated host miRNAs dependent of KSHV infection, indicating that the presence of KSHV has a significant impact on the metabolism of host miRNAs, whose target genes are mainly associated to angiogenesis, ECM, spliceosome, p53 signaling, viral infections, and cell cycle control. Similarly, functional analysis of KSHV miRNA targets showed enrichment in processes, such as cell cycle, spliceosome, RNA transport, microRNA regulation of DDR, and p53 signaling, suggesting that viral miRNAs might mimic cellular miRNAs.
The integrative analysis of viral and host non-coding and coding RNAs and the related processes showed a landscape of the potential relationships of lncRNA-miRNA-mRNA in a KSHV setting. This network highlights that the upregulated genes are involved in processes previously related to KSHV tumorigenesis while downregulated genes are associated with host cell cycle checkpoints and RNA surveillance pathways: G1 to S cycle control, p53 activity regulation, MicroRNA regulation of DDR, Spliceosome, RNA transport, E2F transcription factor network. Finally, the ncRNA-mRNA analysis in the animal model presented here allowed us to develop signatures that can be used to identify druggable gene or networks defining relevant AIDS-KS therapeutic targets.

RNA-Sequencing Analysis
RNA-sequencing raw data used in the present study were obtained as previously described (13). Data are available at https://www.ncbi.nlm.nih.gov/geo/, GSE144101. Briefly, RNA was isolated and purified using the RNeasy mini kit (Qiagen). RNA concentration and integrity were measured on an Agilent 2100 Bioanalyzer (Agilent Technologies). Only RNA samples with RNA integrity values (RIN) over 8.0 were considered for subsequent analysis. mRNA from cell lines and tumor samples were processed for directional mRNA-sequencing library construction using the Preparation Kit according to the manufacturer's protocol. Paired-end sequencing using an Illumina NextSeq500 platform was used, all samples were processed in the same sequencing run of Illumina NextSeq 500 system and analyzed together with the aim to avoid the batches effect. The short sequenced reads were mapped to the mouse reference genome (GRCm38.82) by the splice junction aligner TopHat V2.1.0. Several R/Bioconductor packages to accurately calculate the gene expression abundance at the whole-genome level using the aligned records (BAM files) were used. The number of reads mapped to each gene based on the Mus musculus genome assembly GRCm38 (mm10) were counted, reported and annotated using the featureCounts package. To identify DE genes between cell lines and tumor samples, we utilized the DESeq2 package in R/Bioconductor. DESeq2 performs an internal normalization where geometric mean is calculated for each gene across all samples. The counts for a gene in each sample are then divided by this mean. For ncRNA annotation we employed biomaRt package in R/Bioconductor. We considered the Ensemble transcript ID, the Ensembl gene ID, the Entrezgene ID, the HGNC symbol, the Refseq ncRNA ID and the ReqSeq ncRNA predicted ID. After Deseq2 analysis on all ncRNAs, we filtered out those belonged to the following classes: small nuclear RNA (snRNA), small nucleolar RNA (snRNA), predicted and or pseudogenes, and RIKEN genes; and kept the classes lncRNA and miRNA. originated from frozen batches of mECK36 cells previously generated (12). Briefly, mECs were obtained from Balb/C An Ncr-nu mice (NCI, Bethesda, MD) bone marrow. Mice femurs were flushed twice with PBS, and the eluates were incubated in DMEM media plus 30% FBS (Gemini Bioproducts, Calabasas, CA), 0.2 mg/ml Endothelial Growth Factor (EGF) (Sigma, Saint Louis, MO), 0.2 mg/ml Endothelial Cell Growth Factor Supplement (ECGS) (Sigma, Saint Louis, Missouri), 1.2 mg/1 heparin (Sigma, Saint Louis, MO), insulin transferrin selenium (Invitrogen,Carlsbad, CA), 1% penicillin-streptomycin (Invitrogen, Carlsbad, CA), and BME vitamin (VWR Scientific, Rochester, NY). Cells transfected with KSHVBac36, the vector containing the insert with the genome of KSHV in Bacterial Artificial Chromosome (KSHVBac36) was obtained as described previously, were selected with Hyg-B (12). KSHV (+) tumors were obtained as previously shown, 1x106 KSHV (+) cells were injected subcutaneously into the flanks of nude mice and KSHV (+) tumors formed 5 weeks after injection. KSHV (−) cells were used from frozen populations of KSHV null mECK36 previously obtained (12). KSHV (−) tumor cells were obtained from frozen stocks previously generated by explanted mECK36 tumor cells that have lost the Bac36-KSHV episome (12). These KSHV-negative cells were obtained from frozen stocks previously generated (14). KSHV (−) tumors were obtained as previously shown (12), 1x10 6 KSHV (−) tumor cells were injected subcutaneously into the flanks of nude mice and KSHV (−) tumors formed 3 weeks after injection.

KSHV lncRNA Analysis
Cells and tumors employed in the present study were the same as previously described (13). For lncRNA analysis we included the generated BAM files from eight samples (2 KSHV (+) cells and six KSHV (+) tumors). Based on the KSHV 2.0 reference genome and genome coordinates, we annotated 12 lncRNAs. For measuring gene expression, we applied featureCounts function of the RSubread package in R/Bioconductor. For DEG analysis, we employed DESeq2 package in R/Biocoductor.

Small RNA Sequencing and miRNA Analysis
RNA was isolated and purified using RNeasy Plus Mini Kit (Qiagen, #74134) following the RNeasy MinElute Cleanup Kit (Qiagen, #74204) to separate purification of small RNA (containing miRNA) and larger RNA, the small RNA eluate is enriched in various RNAs of <200 nucleotides. A total of 15 small RNA ranged from cell lines to primary mouse tumors in the presence or absence of KSHV, were processed and sequenced on a HiSeq 2500 System (Illumina, USA). Each sample yielded, on average, 17 million reads, with the exception of one sample (DS016) that was excluded from the analysis for presenting a low number of total reads. Nearly all bases showed scores > Q30 for all reads. Trimmomatic was used to remove adapters and quality control was checked with FastQC. Reads were mapped to a combined mouse and KSHV genome using the bowtie aligner (ver. 1.1.1). To identify novel and known miRNAs we used miRDeep2 package (ver. 2.0.0.7). A hybrid genome of the mouse and the KSHV virus was used for all analyzes in order not to bias the mapping results for or against any of the two separate genomes. The source for all known miRNAs was miRBase (ver. 21). KSHV transcriptome was analyzed using previous resources and KSHV 2.0 reference genome. To identify DE miRNAs across the different comparisons, we utilized the DESeq2 test based on the normalized number of counts mapped to each miRNA. For data integration and visualization of DE transcripts we used R/Bioconductor. Data were submitted to the SRA database, reference PRJNA602753.
Functional enrichment analyses were performed using the ClueGo Cytoscape's plug-in (http://www.cytoscape.org/) and the Enrichr resource (https://maayanlab.cloud/Enrichr/) based on the lists of EVT that were in turn deregulated transcripts across the different comparisons of our model. For pathways terms and annotation, we used those provided by KEGG and BioPlanet (http://tripod.nih.gov/bioplanet/; https://www.genome.jp/kegg/ pathway.html). Significant pathways were based on the Bonferroni Adjusted p value <0.05. To combine and integrate expression data with the results of the functional analysis we used the GOplot package. For the construction of the networks, we used Sankey plots.
All statistical analyses and data visualization plots were done with R/Bioconductor packages.

Genome-Wide Analysis of Non-Coding RNAs in a Cell and Animal Model of Kaposi's Sarcoma
To analyze the ncRNA expression profile in a cell and animal model of Kaposi's sarcoma, we performed deep RNA-seq analysis of all the stages of this model. Figure 1A shows a schematic representation of the model: tumors formed by KSHV Bac36 transfected mouse endothelial cells, KSHV (+) cells, are all episomally infected with KSHV Bac36, and when KSHV (+) cells prior to form tumors lose the KSHV episome in vitro by withdrawal of antibiotic selection, KSHV (−) cells, they completely lose tumorigenicity (12,13). In contrast to KSHV (−) cells, cells explanted from KSHV (+) BAC36 tumors and grown in the absence of antibiotic selection lose the KSHV episome, KSHV (−) tumor cells, are tumorigenic and are able to form KSHV (−) tumors (12)(13)(14).
Unsupervised clustering ( Figure 1B) and Multidimensional plot ( Figure 1C) of the host ncRNAs shows how KSHV status and tissue type cluster with each other. As was previously reported based on mRNA proliles, in vitro and in vivo models clustered separately (13).
To identify changes in host lncRNAs expression profile, we analyzed the number of differentially expressed (DE; FC>1.5, p value <0.05) lncRNAs in key biological comparisons that were detected by RNA-sequencing analysis of: two KSHV (+) cells, two KSHV (−) cells, six KSHV (+) tumors, two KSHV (−) tumor cells and three KSHV (−) tumors (Supplementary Table 1). This mouse model allows for unique experimental comparisons in the same cell and KS-like mouse tumor types: 1) KSHV (−) cells versus KSHV (+) cells can be used to study KSHV mediated effects in vitro, 2) KSHV (−) tumors versus KSHV (+) tumors can be used to dissect the role of ncRNAs in tumorigenesis by comparing tumors driven by KSHV versus tumors driven by host mutations, 3) KSHV (+) cells grown in vitro and in tumors can be used to study in vitro versus in vivo variations induced by micro-environmental cues, and 4) KSHV (−) tumor cells versus KSHV (−) tumors can be used to study in vitro versus in vivo variations in the absence of KSHV (13). We first analyzed lncRNAs expression in these comparisons and found that the highest number of DE lncRNAs was observed in KSHV (+) tumors in both comparisons versus KSHV (−) tumors and versus KSHV (+) cells ( Figure 1D and Supplementary Table 1).

Identification of Relevant lncRNAs in KSHV (+) Tumors
We performed heat map representations of all or top-50 DE lncRNAs -according to each comparison-for all the four biological relevant comparisons mentioned previously (

Pathway Analysis of the lncRNAs, Reveals KSHV Closely Related Bioprocesses
To contextualize the selected lncRNAs into functional processes, we employed LncRNA2Target v2.0 and LncTarD databases (15,16) as resources of lncRNA-target relationships. Since the four selected lncRNAs have been studied more extensively in humans than in mice we searched for their experimentally validated targets (EVT) (Supplementary Table 2). Functional enrichment analysis (KEGG) of the resulting lists of genes revealed several related pathways common to the four lncRNAs, mainly cancer-related pathways and bioprocesses associated with viral diseases (Supplementary Figure 1 and Supplementary Table 2). Interestingly, KSHV infection and MicroRNAs in cancer were the common signature of the 4 lncRNAs (Supplementary Figure 1).
Next, we established a list of the total human target genes contributed by the 4 lncRNAs and looked for their homologues among the DE host genes previously obtained across the different comparisons of our model ( Figure 3A and Supplementary  Table 3). Figures 3B-D, shows the chord plots illustrating the biological process terms and the target genes of the four lncRNAs contributing to that enrichment arranged in order of their expression level in the corresponding comparisons (Supplementary Table 3). Processes such as Integrins in Angiogenesis (with the genes Spp1, Vegfa, Fn1, Kdr, Igf1r), Signaling by PDGF (Vegfa, Kdr, Cdkn1a, Igf1r), HIF-1 signaling (Vegfa, Hif1A, Stat3, Il6, Cdkn1a, Pik3r1, Igf1r), MicroRNAs in cancer (Dicer1, Zeb1, Zeb2, etc.) and KSHV infection (Fgf2, Hif1a, Stat3, Il6, Pik3r1, Jak2, Rb1, Jun) were overrepresented by upregulated target genes in KSHV (+) tumors compared with KSHV (−) tumors. Apoptosis (Mdm2, Bax, Myc, Casp3) and p53 activity regulation (Mdm2, Bax, Casp3, Pcna) were instead associated with downregulated genes in the KSHV-bearing tumors ( Figure 3B and Supplementary Table 3). Similar findings were observed in the comparison KSHV (+) cells  Table 3). Also, pathways of DNA integrity control and cell cycle checkpoints, such as Tp53 network, MicroRNA regulation of DDR and G1 to S cell cycle control were revealed in this comparison, represented by the downregulated genes E2f1, Bax, Dnmt1, Cdkn1a, Myc, or Mdm2. Interestingly, in the transition in vitro to in vivo but in the absence of KSHV we found processes related to the terms KSHV infection (Ccnd1, Ctnnb1, Map2k2, Stat3, Pik3r1, Jun, Erbb2, Igf1r) and MicroRNAs in cancer (Dicer1, Ccnd1, Ctnnb1, Sp1, etc.) associated with downmodulated genes in KSHV (−) tumors ( Figure 3D and Supplementary Table 3), in contrast to that observed in the KSHV-dependent transition ( Figure 3C). Taking together, the integrative in-silico analysis of the lncRNAs-EVT and their associated pathways, with the host transcriptome derived from our model, reveals that the upregulation of Malat1, Neat1, H19, and Meg3 in KSHV (+) tumors would contribute with a network of gene-pathways closely related with KSHV oncogenesis.

KSHV-Dependent In Vitro to In Vivo Transition Is Defined by a Significant Up-Regulation of Host miRNAs
LncRNAs have been demonstrated to regulate gene expression by various mechanisms, including epigenetic modifications, lncRNA-miRNA specific interactions, and lncRNAs as miRNA precursors. Our previous approach showed clear relationship among the four selected lncRNAs and miRNAs in cancer. Therefore, we performed small-RNA sequencing on the samples obtained from our model to identify host DE miRNAs. Next, we conducted an integrated bioinformatics workflow to elucidate relevant networks of lncRNA-miRNA-mRNA during KSHV tumorigenesis.
Unsupervised analysis of 14 samples based on miRNAs expression profiles shows how they cluster together in an unsupervised way according to their predefined features ( Figure 4A). Interestingly, the samples cluster in the same pattern as when the analysis was made for lncRNAs ( Figure 1B) and also for all host genes in our previous study (13). The distance among groups is reflexed in the number of DE miRNAs ( Figure 4B and Supplementary Table 4). Remarkably, the higher proportion of upregulated miRNAs was observed in KSHV (+) tumors (95% of DE miRNAs) compared with KSHV  Table 4). This result is consistent with that previously described in which the term MicroRNAs in cancer was associated with upregulated genes in KSHV (+) tumors and down-regulated genes in KSHV (−) tumors ( Figures 2B, C). Such difference could be partly explained by Dicer1, a master regulator of miRNA biosynthesis, which is in turn linked to the lncRNA H19 ( Figures 3B, D). We performed heat map representations of all or top-50 DE microRNAs -according to each comparison-for all the 4 biological relevant comparisons mentioned previously ( Figure 4C).

Differentially Expressed miRNAs Regulate Gene Targets Related to KSHV Affected Biological Processes
Mature miRNAs regulate gene expression at the posttranscriptional level via partial base-pairing with their target mRNAs. Such interaction leads to mRNA degradation and/or translational inhibition, causing the downregulation of proteins encoded by the miRNA-targeted mRNAs, a biological phenomenon termed RNA interference (RNAi) (21). In silico-based functional analysis of miRNAs usually consists of miRNA target prediction and functional enrichment analysis of miRNA targets.
In the in vitro to in vivo transition, with a more proportional distribution of DE miRNAs, upregulated and downregulated target genes were consequently identified, which provided greater enrichment of bioprocesses closely related to the obtained with the lncRNA targets in the same comparison. As can be seen in the chord plot of Figure 5C, upregulated target   Figure 5C and Supplementary Table 5).
Within the latter, it is worth highlighting the presence of numerous genes related to the nuclear export machinery (Ranbp1, Ran, Xpo1, Nup50, Nup153, Kpnb1, Ncbp1). Lastly, in the in vitro to in vivo transition in the absence of KSHV, fewer terms were significantly over-represented by the miRNAs target genes. Among them highlights Tight junctions, linked to downregulated genes, and Arf6 signal transduction over-represented by the up-regulated ones ( Figure 5D and Supplementary  Table 5).
Collectively, these results indicate that the presence of KSHV has a significant impact on the metabolism of host miRNAs, which contribute to the regulation of host genes linked to processes of angiogenesis, ECM, transcriptional metabolism, viral infections and cell cycle control, mainly.
Among other non-coding RNAs KSHV encodes a number of lncRNAs (7). We inquired into the RNA-seq data and identified seven out of twelve annotated lncRNAs, with detectable levels of expression in KSHV (+) cells and tumors. As-ORF7, as-K5/K6, as-ORF65/69, and ALT were the most abundant transcripts (Supplementary Table 6

Identification of a lncRNA-miRNA-mRNA Interaction Network Involved With KSHV Tumorigenesis
LncRNAs can also serve as regulatory elements of the RNAi pathway (22). Indeed, host lncRNA transcripts are involved not only with the maturation of miRNA transcripts but also they may interfere with miRNA induced translation inhibition, thus acting as competing endogenous RNAs (ceRNAs), or "sponge RNAs" (22). Such lncRNA-miRNA associations allow for a fine tuning of gene expression regulation. Therefore, dysregulation of the lncRNA-miRNA balance could contribute to the onset of KSHV pathogenesis.
To identify relevant pairs of lncRNA-miRNA in our model, we used DIANA-LncBase v3 (18). For each of the four lncRNAs we searched for their highly confident experimentally supported viral and host miRNA targets, derived from highthroughput methodologies, which were in turn DE in the corresponding comparison.
Within the ten most abundant KSHV miRNAs, we found that five of them have been associated with the human lncRNAs MALAT1 and NEAT1. When analyzing the targets (downregulated in KSHV + tumors) of these five miRNAs, we observed that they share most of the genes obtained with the ten miRNAs, which is therefore reflected in the same enriched pathways ( Figure 6C and Supplementary Table 7). Using the same approach, we proceeded with the analyses for the lncRNA-host miRNAs associations. For the in vitro to in vivo transition dependent of KSHV, 31 miRNAs accomplished the criteria, of which 23 were upregulated and 8 were downregulated in KSHV (+) tumors (Figure 7 and Supplementary Table 7). The highest contribution was made by Malat1 (29 miRNAs) and Neat1 (21 miRNAs), followed by Meg3 (11 miRNAs) and H19 (8 miRNAs). Among the downregulated miRNAs highlights the members of the miR17-92 family: miR-17-5p, miR-19a-3p, miR-20a-5p, and miR-92a-3p. Their respective targets are represented by genes such as Egfr, Foxo1, Pdgfra, Rb1, Igf1, Map3k1, etc. all upregulated in KSHV (+) tumors (Figure 7). Other relevant downregulated miRNAs were miR-128-3p and miR-155-5p, which target multiple common genes. On the other hand, upmodulated miRNAs were linked mostly to Malat1 and Neat1.
Remarkably, among them are miR27-b-3p, miR-140-3p, miR-142-3p, and miR-142-5p, whose gene precursors were also found up-modulated in KSHV (+) tumors ( Figure 2B). As can be seen in Figure 7, the functional analysis that arose from the lncRNA-miRNA-mRNA triad shows that the pathways are arranged in an unsupervised way in three main clusters. The MAPK signaling together with Pathways in cancer would make up the 1st group, over-represented by the upregulated target genes contributed mainly by miR-671-5p, miR-128a-3p, miR-155-5p, and let 7e-5p, as well as the miRNAs of the miR17-92 cluster. A second group is integrated by processes related to Viral infection (HPV infection), Matrix organization and Angiogenesis, represented by genes contributed by the miRNAs of the miR17-92 cluster, let-7d-5p and miR-124-3p, along with others ( Figure 7). A third group would be made up of the pathways HIV life cycle and p53 signaling, over-represented by negatively regulated genes, targets of the miRNAs miR-27b-3p, miR-101a-3p, miR-140-3p, and miR-142, among others (Figure 7).
For the comparison KSHV (−) tumors versus KSHV (+) tumors, we obtained a network of the four lncRNAs targeting 26 miRNAs all upregulated in KSHV (+) tumors with their corresponding downregulated target genes (Figure 8 and Supplementary Table 7). It is evident a shift in the expression of specific miRNAs, such as let-7e-5p, let-7d-5p, miR-123-3p, and miR-31-5p, compared to that observed in the KSHV-dependent transition. Other relevant miRNAs that appear are miR-26b-5p, miR-181 (with its variants a, b and c), miR-378-3p, and miR-381-3p. Here again, the presence of miR-140-3p and miR-378-3p correlates with their respective immature precursors that had been identified previously as upregulated along with the lncRNAs ( Figure 3B). By analyzing the mRNA targets of the miRNA signature, previously identified as downregulated in KSHV (+) tumors, we obtained a relatively small group of genes that function in two major related processes: the regulation of cell cycle control (G1 to S cycle control, p53 activity regulation, MicroRNA regulation of DDR) and the transcription machinery, with the pre-mRNA splicing machinery (Spliceosome) and the E2F transcription factor network ( Figure 8). Remarkably, this functional pattern resembles that observed with the KSHV miRNAs ( Figure 6).
Collectively, our analysis revealed a functional network of lncRNA-miRNA-mRNA in a KSHV animal model.

Gene Signatures Used to Identify Drug-Associated Genes or Networks
KS remains potentially life threatening for patients with advanced or ART-resistant disease, where systemic therapy is indicated and three FDA-approved agents that include liposomal anthracyclines are available (4,23,24). Despite the effectiveness of these agents, most patients progress within six to seven months of treatment and require additional therapy (25). Therefore, there is a need to develop alternative strategies. Identifying drugs or clinical candidates that synergize with the current KS frontline therapeutic approaches has immediate translational potential that would be realized in a clinical trial if identified drug combinations show sustained efficacy in animal models. Our animal model allowed us to develop signatures that can be used to identify druggable gene or networks defining relevant AIDS-KS therapeutic targets.
For this end we used two approaches: 1) druggable miRNAsgene pairs, and 2) the complete signature of the lncRNA-miRNA-mRNA network for upregulated genes.
Importantly, some of the aforementioned drugs (abacavir, doxorubicin, bevacizumab, bortezomib, imatinib, sirolimus, and thalidamide) have been evaluated alone or in combination with other drugs in different KS clinical trials. The description of such studies is found in Supplementary Table 8. The fact that our analyses pointed to drugs that target KS oncogenic pathways identified in the laboratory or drugs that are currently in use of being tested in AIDS-KS, reinforces the possibility of involvement of the KSHV regulated ncRNA network in viral sarcomagenesis.

DISCUSSION
Virus-host interactions trigger a set of mechanisms that eventually affect the expression of host genes involved in the regulation of the viral replicative cycle as well as the pathogenesis of the disease (27). Whereas dysregulation of host protein-coding genes caused by KSHV infection is well explored, host ncRNAs and KSHV dependency remains poorly characterized. Currently, miRNAs and lncRNAs are by far two of the most commonly studied ncRNA biotypes (28,29). We have previously developed and characterized a unique multistep KSHV tumorigenesis model in which cells explanted from a KSHV (+) tumor that lose the episome can form KSHV (−) tumors driven by host mutations such as the PDGFRA-D842V (12,30). Using NGS on this model, we interrogated the transcriptional, genetic and epigenetic (CpG island methylation) landscape upon KSHV tumor formation and upon KSHV-loss in cells and tumors (13). In such study, we focused on the host and virus coding genes. Therefore, taking advantage of the model and the RNA-sequencing technology, we decidedfor this study-to explore the transcriptional consequences of KSHV tumorigenesis on the ncRNAs setting, with the aim of identifying a functional interplay between lncRNAs and miRNAs dependent of KSHV.
Here we identified four relevant lncRNAs upregulated in KSHV (+) tumors: Malat1, Neat1, H19 and Meg3. Accumulating evidence has shown that lncRNA exert its functions by regulating the expression of target genes. As a first approach, using databases that collects all lncRNA-target relationships confirmed by binding experimental technologies, we searched for the target genes for the human homologues of each of the selected lncRNAs. In addition to having common target genes, pathway analysis showed that the four lncRNAs also share common related processes, mainly associated with cancer and viral infections. Interestingly, KSHV infection and MicroRNAs in cancer were among the common overrepresented terms.
Next, we interrogated the transcriptome of our model to identify the 4-lncRNAs common targets into the DEG. The integrated analysis allowed us to define a reduced group of host lncRNAs-target genes that significantly would contribute with KSHV tumorigenesis and related processes. The integration of the in silico approach of the lncRNAs-EVT and their associated pathways, with the host transcriptome derived from our model, reveals a network of gene-pathways closely related with KSHV oncogenesis: Integrins in angiogenesis, KSHV infection, signaling by PDGF, HIF1-signaling pathway or MicroRNAs in cancer were represented by upregulated genes such as Egfr; Vegfa, Hif1a, Dicer1, Zeb1, Zeb2, Rb1, or Il6.
In addition, one of the distinctive pathways of the in vitro to in vivo transition dependent of KSHV, provided by the lncRNA targets, was Extracellular Matrix Organization and Activation of Matrix Metalloproteinases, overrepresented by the MMPs Mmp2, Mmp9, Mmp13, and Mmmp14. MMPs are associated with KS and may contribute to the mechanism of KS tumor growth. They are usually synthesized by the tumor stromal cells, including fibroblasts, myofibroblasts, inflammatory cells and endothelial cells. These components can also integrate a tumor derived from cells in vitro. Although the mechanism by which Malat1, Neat1, or H19 regulate the expression of MMPs is not yet clear, different studies have shown that the silencing or overexpression of these lncRNAs positively correlate with the expression of MMPs, such as MMP9 or MMP2 (31)(32)(33).
MALAT1 is perhaps the most studied lncRNA and consequently the one with the most targets. It has been shown to regulate EGFR expression promoting carcinogenesis (34); it has been shown to regulate endothelial cell function and vessel growth (35); it has been defined as a hypoxia-induced lncRNA (36); it modulates ZEB1 and ZEB2 by sponging miRNAs (37,38). Remarkably, MALAT1 expression is induced by the plateletderived growth factor BB (PDGF-BB) (39). In a recent study, we have shown that the KSHV-ligand mediated activation of the PDGF signaling pathway is critical for KS development (30). Later, we found that two PDGFs, Pdgfa and Pdgfb, and their receptor Pdgfra were both hypo-methylated and up-regulated in KSHV (+) tumors (13). Overall, the evidence clearly shows that Malat1 is a key regulator of several target genes involved in KSHV-dependent signaling pathways. It remains to be determined whether Malat1 is a driver or simply a passenger of KSHV tumorigenesis.
NEAT1, is closely related to MALAT1 (aka NEAT2), and both have been shown to bind multiple genomic loci on active genes, but display distinct binding patterns, suggesting independent but complementary functions (40). As MALAT1, NEAT1 is retained in the nucleus where it forms the core structural component of the paraspeckle sub-organelles. The formation of paraspeckle increases in response to viral infection or proinflammatory stimuli (41). Furthermore, Viollet et al. (42) demonstrated that NEAT1 is upregulated in KSHV infected cells versus noninfected cells under hypoxic conditions. Our results show that Neat1 is upregulated in KSHV-cells versus KSHV+ cells and indeed is upregulated in KHSV (+) tumors, during the in vitro to in vivo transition. On the other hand, the lncRNA target analysis showed that Neat1 positively associates with the upregulated targets Il6, Stat3 and Spp1 in the KSHV (+) tumors. In this regard, NEAT1 has been shown to strengthen IL-6/STAT3 signaling and promote tumor growth and proliferation through nuclear trapping of mRNAs and proteins which acts as inhibitors of the IL-6/STAT3 signaling pathway (43). Previously, it had been demonstrated that STAT3 is activated by KSHV infection and correlates with IL6 release in dendritic cells (44). In summary, these data taking together reveal a host network in which upregulation of Neat1 would favor the activation of IL6/STAT3 signaling contributing directly or indirectly to KSHV tumorigenesis. MEG3 is generally considered as a tumor suppressor lncRNA. In this study we found a downregulation of Meg3 in KSHV (−) cells versus KSHV (+) cells. However, a significant increase of the lncRNA was evidenced in the in vitro to in vivo transition. Sethuraman et al. (11) showed that KSHV employs its miRNAs to target MEG3 promoting its downmodulation to potentially contribute to sarcomagenesis. Therefore, it is possible to speculate on a downmodulation of Meg3 by the expressed KSHV miRNAs as an early event in the viral cycle followed by an upmodulation of Meg3 as a response of the host cell to the already triggered tumor growth.
KSHV drives latently infected cells towards proliferation by a variety of mechanisms such as interfering with MEG3 or the p53 pathway through miRNAs or the protein LANA, respectively (11,45). In this study, we identified that p53 network would be regulated in a KSHV-dependent manner by the modulation of key genes targeted by the lncRNAs, such as Casp3, Bax, Mdm2, Cdkn1a, or Pcna. Interestingly, these genes along with E2f1 are linked to other related processes such as G1 to S phase regulation and MicroRNA DDR. KSHV needs to face various cellular defense mechanisms designed to eradicate the viral infection. One such response can include DDR response factors, which can promote an arrest in cell growth (G1-S regulation) and trigger cell death (p53 network, Apoptosis). Our findings indicate that those processes would be repressed through the downmodulation of the mention lncRNA targets in KSHV (−) tumors versus KSHV (+) tumors, as well as in the KSHV in vitro to in vivo transition. Remarkably, several studies have shown that viruses including KSHV have developed suppressive strategies against DDR (9,46). In this sense, KSHV miRNAs are relevant for protecting cells from DDR (47,48). In addition, cellular lncRNAs are important gene regulators of DDR in a process which involve essential players of miRNA biosynthesis such as DICER1 and DROSHA (22,48). In fact, Dicer1 was one of the significantly upregulated target genes linked to H19 in the comparison KSHV (−) tumors versus KSHV (+) tumors. In summary, there is a complex network between KSHV and host ncRNAs that would regulate DDR factors in order to bypass cell cycle checkpoints. miRNA analysis revealed a high proportion of upregulated host miRNAs dependent of KSHV infection. This finding led us to interrogate the functional processes associated to the miRNAs targets. Enriched terms were linked to p53 signaling, Spliceosome and Cell cycle. When evaluating the in vitro to in vivo transition which involved both up and downregulated miRNAs, processes such as Integrins in Angiogenesis, Platelet activation, or signaling by PDGF were associated to the upregulated targets, whereas viral infection (HPV, HIV) or p53 signaling were linked to the downregulated targets.
Furthermore, we evaluated the relevance of viral lncRNAs and miRNAs expression in KSHV tumorigenesis. KSHV encodes at least 16 potential lncRNAs (49). In our analysis, we were able to annotate 12 lncRNAs of which 7 showed detectable levels of expression in KSHV (+) cells and tumors. PAN RNA (polyadenylated nuclear RNA), the most abundant and characterized KSHV lncRNA linked to KSHV lytic gene expression (49), is expressed in KSHV (+) tumors (Supplementary Figure 2) correlating with the in vivo up-regulation of KSHV lytic gene expression (13). In addition, we identified as-ORF7 and as-K5/6 upregulated in KHSV (+) tumors compared to KSHV (+) cells. Although their functions are still not reported, our results indicate that these transcripts would have a potential role in KSHV tumorigenesis. Regarding miRNAs, we identified a group of ten relevant members which constituted the most frequent in mouse KSHV (+) tumors. Among them highlights K12-4-3p, K12-3-5p, K12-8-3p, previously identified as highly expressed in human KS lesions (50). Moreover, K12-4-3p, which represented 50% of the KSHV miRNAs detected in this analysis in mouse KSHV (+) tumors, was shown to be able to restore the transforming phenotype of a mutant KSHV containing a deletion of all KSHV microRNAs (51), indicating its association with cellular transformation and tumor induction. Similarly, K12-3-5p was shown to promotes cell migration and invasion of endothelial cells (52). The functional analysis of their targets-downregulated in mouse KSHV (+) tumors-showed enrichment in processes such as Cell cycle, Spliceosome, RNA transport, MicroRNA Regulation of DDR, and p53 signaling, coinciding with what was observed with the host miRNAs, which suggests that viral miRNAs might mimic cellular miRNAs. It is possible that the same targets are also relevant to the infection of human cells by KSHV (KSHV miRNA) and to KSHV pathogenesis (host miRNA) (53,54).
Since the similarity with the one found with respect to the lncRNA targets, we decided to carry out an integration network dependent of KSHV between the 4-lncRNAs, the DE miRNAs (from virus and host) related to those lncRNAs, their validated targets, and the related processes.
The integration showed a more concise landscape of the potential relationships of lncRNA-miRNA-mRNA in a KSHV setting, in which, once more, highlights that the upregulated genes are involved in processes, such as pathways in cancer and those previously closely related to KSHV tumorigenesis, including Angiogenesis, PDGF signaling, MAPK signaling or ECM organization. Similarly, down-modulated genes linked preferentially to p53 signaling, Spliceosome, miRNA regulation of DDR, or RNA transport, among others. In the latter Proteasome subunit protein PSMD1, nucleoporins NUP50, NUP153, nucleolar protein NPM1, or exportin 1 XPO1 have been shown to modulate HIV infection or other viral cycles (55)(56)(57)(58). In addition, in a preprint article it was postulated an extensive destruction of the nuclear and nucleolar architecture during lytic reactivation of KSHV, with redistribution or degradation of proteins such as NPM1 (59). More interestingly NPM1 (aka NPM) is a critical regulator of KSHV latency via functional interactions with v-cyclin and LANA. Strikingly, depletion of NPM in PEL cells has led to viral reactivation, and production of new infectious virus particles (60). On the other hand, using a model of oncogenic virus KSHV-driven cellular transformation of primary cells, Gruffaz et al. (61) illustrate that XPO1 is a vulnerable target of cancer cells and reveal a novel mechanism for blocking cancer cell proliferation by XPO1 inhibition.
Spliceosome has been other of the relevant terms yielded by our network analysis. In the presence of KSHV, positively regulated miRNAs linked to a group of down-modulated targets closely related to the splicing machinery (Snrpb, Snrpb2, Rbmx, Srsf3, Srsf10). NEAT1 and MALAT1 were the first lncRNAs to be identified as having a relevant role in mRNA splicing in both human and mouse cells. The mechanisms by which both lncRNAs modulate splicing is extensively reviewed in Romero-Barrios et al. (62). Remarkably, it has been postulated that MALAT1 modulates the phosphorylation status of a pool of Serine/arginine-rich (SR) proteins (proteins involved in splicing), resulting in the mislocalization of speckle components and changes in alternative splicing of pre-mRNAs, impacting in other SR-dependent post-transcriptional regulatory mechanisms, including RNA export, NMD and translation (63). In addition, cells depleted for MALAT1 show an increased cytoplasmic pool of poly(A)+ RNA, suggesting that MALAT1 contribute with the retention of nuclear mRNAs. Before RNAs can interact with nuclear export machinery, they must undergo processes that regulate the number of transcripts that is exported to the cytoplasm or nuclear decay pathways. KSHV manipulation of nuclear RNA regulation is one of the strategies acquired by the virus to influence the host RNAs during viral infection (8). In fact, it was very recently demonstrated that NMD pathway targets KSHV RNAs to restrict the virus (10). In summary, our network reveals another intricate relationship between lncRNA-miRNAtargets that can function in modulating spliceosome pathway and RNA transport during virus-host interaction.
Other relevant miRNAs that emerged from our network were members of the cluster 17-92 and the let-7 family whose multiple targets regulate different pathways associated with cancer. Moreover, miR 140-3p or miR378b also stand out, of which, like miR143-3p, their precursors were found upregulated in KSHV (+) tumors.
It has been demonstrated that miR140 in the nucleus can interact with NEAT1, leading to the increased NEAT1 expression (64). Remarkably, there is another interesting link with NEAT1, which is p53 signaling, a frequent pathway represented in our networks. It has been shown that silencing Neat1 in mice prevents paraspeckle formation, which sensitizes preneoplastic cells to DDR activating cell death and impairing skin tumorigenesis (65). Moreover, activation of p53 stimulates the formation of NEAT1 paraspeckles, establishing a direct functional link between p53 and paraspeckle biology (65). P53 regulates NEAT1 expression to stimulate paraspeckle formation and NEAT1 paraspeckles, in turn, dampen replication-associated DNA damage and p53 activation in a negative regulatory feedback (65). These data indicate that upregulation of Neat1 in KSHV (+) tumors could attenuate p53 signaling network, and infected cells may benefit from this situation evading the p53 checkpoint in response to DNA damage.
As mentioned before, we have used this same mECK36 tumor model to analyze the consequences of KSHV loss by comparing the mutational and methylation landscape of KSHV (+) and KSHV (−) tumors. We found that KSHV loss led to irreversible oncogenic alterations including oncogenic mutations and irreversible epigenetic alterations that were essential in driving oncogenesis in the absence of KSHV (13). In contrast to these irreversible effects of KSHV tumorigenesis, the ncRNA network we describe in the present study display a high degree of plasticity and reversibility upon KSHV loss further supporting the idea that these oncogenic networks are driving tumorigenesis and are more strictly dependent on the presence of KSHV.
Along the study, we integrate mice genes and their homologues in humans to understand the ncRNA biology in KSHV tumorigenesis and to develop signatures that can be used to identify druggable gene or networks defining relevant AIDS-KS therapeutic targets.
Interestingly, we identified drugs usually used against targets in experimental KSHV models or in clinical trials: Abacavir, Bevacizumab, Bortezomib, Celecoxib, Doxorubicin, Imatinib, Oxaliplatin, Sirolimus, Sunitinib, Thalidomide and Vorinostat. Interestingly, we have previously shown a combinatory effect between Bortezomib and Vorinostat for the treatment for primary effusion lymphoma (66). The fact that our analyses pointed to drugs that target KS oncogenic pathways identified in the laboratory or drugs that are currently in use of being tested in AIDS-KS further validate the bioinformatic analysis in our KSHV mouse tumorigenic model and reinforces the idea of the involvement of the KSHV regulated ncRNA network in viral sarcomagenesis.
In summary, in the present study the integration of the transcriptional analysis of ncRNAs in a KSHV model in cells and mouse tumors, with an exhaustive computational analysis of their experimentally supported targets, has allowed us to dissect a complex network that defines the main pathways involved in KSHV pathogenesis and host response. Understanding the relationships between these different RNA species will allow a better understanding of the biology of KSHV and can aid in the identification of relevant AIDS-KS druggable targets.

DATA AVAILABILITY STATEMENT
The data sets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
JN: conceptualization, investigation, methodology, resources, supervision, writing-original draft, writing-review and editing. MS: formal analysis, methodology, writing-original draft. DS: methodology. SR: methodology. SW: methodology. OC: funding acquisition, writing-review and editing. MA: funding acquisition writing-review and editing. EM: funding acquisition, resources, supervision, writing-review and editing. EL: conceptualization, formal analysis, funding acquisition, investigation, resources, supervision, writing-original draft, Writing-review and editing. All authors contributed to the article and approved the submitted version.