Differential expression analysis comparing EBV uninfected to infected human cell lines identifies induced non-micro small non-coding RNAs.

Epstein-Barr virus (EBV) is a ubiquitous human herpes virus, which is implicated in cancer and various autoimmune diseases. This study profiles non-micro small non-coding RNA expression changes induced by latent EBV infection. Using small RNA-Seq, 346 non-micro small RNAs were identified as being significantly differentially expressed between EBV(+) BJAB-B1 and EBV(-) BJAB cell lines. Select small RNA expression changes were experimentally validated in the BJAB-B1 cell line as well as the EBV-infected Raji and Jijoye cell lines. This latter analysis recapitulated the previously identified induction of vault RNA1, while also finding novel evidence for the deregulation of several tRNAs and a snoRNA.

Introduction

Epstein–Barr virus (EBV), also known as human herpesvirus 4 (HHV-4), is a widespread γ-herpes virus that infects over 90% of the human population [1]. EBV, the first ever virus to be associated with cancer, was discovered in 1964 [2] when it was isolated from cultured lymphoblasts derived from Burkitt's lymphoma (BL) patients [3]. Since then, EBV has been associated with various cancers [4,5] as well as autoimmune diseases [6,7].

EBV is estimated to account for more than 200,000 cases of cancer each year and EBV-associated malignancies are responsible for 1.8% of cancer-related deaths [8]. Initial EBV replication, due to infection, occurs via an aggressive lytic phase [9,10]. After replication, the virions enter and reside in long-lived memory B cells causing persistent life-long latent infections with continuous lytic reactivation to produce infectious viral particles [4,11]. EBV infection is associated with a number of malignancies, including a variety of B-cell lymphomas, NK/T-cell lymphoma, nasopharyngeal carcinoma, and gastric carcinoma [4]. In B cells, the best characterized malignancy-associated cell type, EBV progresses through multiple distinct latency states in which small non-coding RNAs (EBERs and miRNAs) are consistently expressed and viral antigen genes are progressively downregulated [12].

Advances in next generation sequencing have contributed to an explosion of knowledge regarding the EBV genome and transcriptome, as well as the effects of EBV on the human genome: discovery of novel small non-coding RNAs expressed by EBV [9], EBV encoded circular RNAs (circRNAs) [13,14], EBV and host miRNA expression in patient and cultured cell samples [[15], [16], [17], [18], [19], [20], [21], [22], [23], [24]], differential expression of host and EBV mRNAs, long non-coding RNAs (lncRNAs) as well as sequence variations in EBV genomes [9,[25], [26], [27], [28], [29]]. However, the effect of infection on host non-micro small non-coding (nc)RNA has not been extensively studied.

In this study, we focused on identifying differentially expressed small ncRNAs between two widely-used laboratory cell lines: EBV(+) BJAB-B1 and EBV(−) BJAB cell lines. BJAB cells were derived from a human BL [30]; subsequent in vitro infection of BJAB cells using the P3HR1 (Type II) strain of EBV was used to generate BJAB-B1 cells [31]. This created a paired set of cells that are isogenic but for the presence of EBV, which facilitates the study of EBV effects on host cells. Additionally, BJAB-B1 cells express a viral latency program (latency III) that recapitulates early points in infection. An analysis of differentially expressed small ncRNAs between these common cell lines will increase understanding of their unique biologies and may contribute to a better understanding of EBV-associated maintenance or progression of malignancy (e.g. via the deregulation of host non-micro small ncRNAs).

Materials and methods

Cell culture and RNA extraction

BJAB, BJAB-B1, Raji, and Jijoye cell lines were grown in RPMI 1640 media, supplemented with 2 mM l-glutamine, 10% fetal bovine serum (FBS), and 1% penicillin-streptomycin, at 37 °C in a 5% CO2 environment. Cells were spun down at 150 x g for 5 min and the resulting pellet was dissolved in 1 mL of cold TRIzol (Invitrogen) and allowed to incubate at room temperature for 5 min. The cell lysate and TRIzol mixture was added to Phase Lock Gel (PLG) tubes (QuantaBio) with 200 μL of chloroform, shaken vigorously for 15 s and incubated at room temperature for 3 min prior to centrifugation at 12000 x g at 4 °C for 5 min. The top aqueous layer was transferred to a new tube and 1 μL of Glycoblue (Invitrogen) and 500 μL of isopropanol was added and incubated at room temperature for 10 min. The RNA was pelleted by centrifugation at 12000 x g for 10 min at 4 °C. The pellet was washed with 1 mL of 70% ethanol, vortexed briefly, and centrifuged at 7500 x g at 4 °C for 5 min. The RNA pellet was dried briefly and resuspended in nuclease free water. This RNA was used for the experimental procedures described below.

Library construction and small RNA sequencing

For small RNA sequencing of BJAB and BJAB-B1 cells, total RNA samples were further purified using the miRNeasy Mini Kit (Qiagen) according to manufacturer's protocol (did not enrich for miRNAs). The total RNA quantity and purity were analyzed using a Bioanalyzer 2100 (Agilent, CA, USA) to confirm RIN numbers were >9.0 before 3 biological replicates per BJAB and BJAB-B1 cell line were sent for library preparation and small RNA sequencing at LC Sciences (Houston, USA). Small RNA enrichment was done by excision of the 15 to 200 nt fraction from a polyacrylamide gel. Approximately 1 μg of enriched RNA was used to prepare a small RNA library according to the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, USA) protocol. Paired-end sequencing (150 bp) on an Illumina Hiseq 2500 was carried out following the vendor's recommended protocol.

miRNA bioinformatics analysis

Raw reads were subjected to an LC Sciences in-house program, ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove adapter dimers, low quality sequences, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. Subsequently, unique sequences with length in the 18–26 nt range were mapped to specific species precursors in miRBase 21.0 by BLAST search to identify known miRNAs and novel 3p- and 5p-derived miRNAs. Length variation at both 3′ and 5′ ends and one mismatch inside of the sequence were allowed in the alignment. The unique sequences mapping to specific species mature miRNAs in hairpin arms were identified as known miRNAs.

Non-micro small RNA bioinformatics analysis

The quality of the paired-end, raw sequences as well as the adapters used was checked by FastQC [32]. Illumina TruSeq small RNA 5′ and 3′ adapters were trimmed using Trim Galore [33] that utilizes Cutadapt [34] to trim adapters. The trimmed sequences were reexamined for quality and trimming efficiency with FastQC. The ‘new Tuxedo’ package comprised of HISAT2, StringTie, and Ballgown was used to align reads to the human genome, assemble non-micro small RNA transcripts, and compute the abundance of transcripts and differential expression (DE) of these transcripts between BJAB-B1 and BJAB cell lines according to the stated protocol [35]. Reads from each sample were mapped to the reference human genome obtained from Ensemble FTP server (ftp://ftp.ensembl.org/pub/release-94/fasta/homo_sapiens/dna/): Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz [36] using HISAT2. Transcripts for each sample were assembled with StringTie [37] (Release 28: GRCh38.p12). Transcript abundances were estimated, and table counts were created for DE analysis between BJAB-B1 and BJAB samples via Ballgown. DE analysis was carried out using R packages Ballgown and genefilter; low-abundance genes were filtered out by removing all transcripts with a variance across samples less than one. After this, transcripts that showed statistically significant differences between BJAB-B1 and BJAB were identified and transcript IDs, taken from RNACentral [38,39] annotations, were appended to the results (Table A. 2).

Quantitative reverse transcription PCR

RNA was extracted according to the manufacturer's instructions using Trizol. Total RNA was DNase (M0303S; NEB) treated and recovered using phenol/chloroform extraction and ethanol precipitation. Reverse transcription (RT) reactions were performed with Superscript III using random hexamer primers (Invitrogen). For the RT reaction, the RNA was diluted tenfold to prepare the 1% input. QPCR was performed using PowerUP Sybr mix and the QuantStudio3 instrument (Thermo Fisher). Gene expression levels were assessed by qPCR using the threshold cycle (ΔΔCT) method; no-RT and no template reactions were included. Primers (forward/reverse) used were as follows: tRNA-Asn: TTCTGTAGCGCGATCGGT/GGAACGAACCACCAACCTTT; tRNA-Gln: TAATGGTGAGCACCCTGGA/AGGTACCACCGAGATTTGAAC; tRNA-His: GCCGTGATCGTATAGTGGTTAG/GATTCGAACCGAGGTTGCT; tRNA-iMet: AGAGTGGCGCAGCGGAA/GGATGGTTTAGATCCATCGACCTCTG; tRNA-Leu: CCGGGATGGCTGAGTGGTTAAG/CACGCGGACACCTGTCCATT; tRNA-Thr: GAGTGGTGCAGCGGAAG/AGGATGGCTTCGATCCATTG; SNORD89: CAGTGTCTCCATCAGCAGTTT/CAGACTAGTGGTTCGCTTCAG; vtRNA1: GGCTTTAGCTCAGCGGTTA/GTCTCGAACAACCCAGACAG.

Data availability

All sequencing data have been deposited in the NCBI Gene Expression Omnibus (GEO): Accession # GSE139425. All other data pertaining to this study are presented below or in the Appendix.

Results

Differential expression of human miRNAs

Although the focus of this study was on the discovery of differentially expressed non-micro small ncRNAs, data from the miRNA fraction were also analyzed (Table A. 1). A total of 100 miRNAs were found to be significantly (p-value < 0.05) differentially expressed between BJAB-B1 and BJAB cells. The majority of differentially expressed miRNAs showed low fold-change (less than 2) or were not highly expressed. Limiting results to highly-expressed miRNAs, where (in either cell type) average reads per sample were >900, and that had fold changes > 2, led to 14 miRNAs (Table 1). Each of these miRNAs were previously reported as being deregulated by EBV infection or EBV-associated pathologies [20,[40], [41], [42], [43], [44], [45], [46], [47], [48], [49]]. Our results are consistent with miRNAs expected to be altered in BJAB-B1 cells. For example, the two highest abundance miRNAs in BJAB-B1 cells are miR-21-5p and miR-146a-5p, which are increased by a factor of ~2-fold and ~16-fold, respectively vs. uninfected BJAB cells. The most down-regulated miRNA in BJAB-B1 cells (vs. BJAB) is miR-200c-3p, which has a fold-change of ~18.

Table 1

Highly expressed miRNAs called as being differentially expressed: fold-change (FC), mean expression data, and p-values are presented.

miRNA	FC	log2(FC)	B1(mean)	BJAB(mean)	p-value
miR-155-5p	37.26	5.22	2831	76	2.04E-04
miR-146a-5p	15.97	4.00	40112	2512	2.32E-03
miR-99a-5p	9.63	3.27	13661	1419	8.08E-03
miR-125b-5p	6.40	2.68	1890	295	7.49E-04
miR-26a-5p	4.38	2.13	8968	2048	2.92E-03
miR-21-3p	3.09	1.63	4101	1329	8.42E-04
miR-192-5p	2.53	1.34	2379	941	6.23E-04
miR-22-3p	2.37	1.25	951	401	1.15E-03
miR-29a-3p	2.37	1.24	2138	903	4.28E-03
miR-21-5p	2.06	1.04	52437	25498	2.94E-03
miR-7-5p	0.34	−1.54	1447	4210	1.85E-02
miR-185-5p	0.28	−1.86	369	1341	3.04E-05
miR-28-3p	0.13	−3.00	138	1101	7.12E-05
miR-200c-3p	0.05	−4.20	115	2120	3.80E-02

Differential expression of human non-micro small ncRNAs

Comparing BJAB-B1 and BJAB cells, 346 non-micro small ncRNAs were called as being significantly (p < 0.05) differentially expressed (Table A. 2). These can be binned into 12 families of RNAs, where the three most populated classes are tRNA, pre-miRNA, and snoRNA (accounting for 73% of the identified RNAs; Fig. 1). The tRNA class comprised the largest group (33%). Increases in non-micro small ncRNAs accounted for ~65% of the results (Table A. 2)— 82% of tRNAs were identified as being up-regulated and 81% of snoRNAs were also identified as being significantly up-regulated.

Fig. 1

Pie chart showing the percentage of 12 classes of non-micro small ncRNAs amongst those called as being significantly differentially expressed (from results in Table A. 2).

Eight of the most highly up-regulated small ncRNAs were selected for additional validation (via RT-qPCR): six tRNAs, C/D Box Small Nucleolar RNA (SNORD)89, and vault (vt)RNA1 (Fig. 2 and Table A. 3). Differential expression was measured comparing uninfected BJAB cells to infected BJAB-B1 cells, as well as two additional EBV-infected cell lines: Raji, and Jijoye. The most up-regulated small ncRNA across all cell types was vtRNA1, which was previously noted in EBV-infected B cells using subtractive hybridization [50]. The next most consistently up-regulated ncRNA was SNORD89. It was the second most up-regulated ncRNA in Raji and Jijoye; however, although highly up-regulated in BJAB-B1 (~45X fold-change), it was on-par with the most up-regulated tRNAs in this cell line.

Fig. 2

RT-qPCR analysis of select small ncRNAs. The relative expression of select small ncRNAs (called as highly up-regulated by RNA-Seq) is shown, comparing BJAB-B1, Raji, and Jijoye samples to BJAB samples (data in Table A. 3). The error is represented as the standard error of 3X biological replicates.

In BJAB-B1 cells, tRNA-Asn (GTT anticodon), tRNA-Gln (CTG anticodon), tRNA-Leu (CAA anticodon), and tRNA-Thr (CGT anticodon) were all consistently up-regulated (>45X fold-change) with tRNA-Asn being highest (~147X fold-change; Fig. 2 and Table A. 3). In Raji cells tRNA-Asn was also highly up-regulated, with fold-change ~60X, whereas tRNA-Gln, tRNA-Leu, and tRNA-Thr were not as highly up-regulated. Although RNA-Seq results found tRNA-His and tRNA-iMet to be up-regulated, this could not be recapitulated in the RT-qPCR data (Fig. 2 and Table A. 3). tRNA-His only slightly changed across samples and, in Jijoye cells, was found to be reduced by 1.22-fold. tRNA-iMet levels were slightly reduced in all cell types, with the highest fold-change being 3.25-fold in Jijoye cells.

Discussion

Small RNA-Seq results for human miRNAs differentially expressed by EBV infection were consistent with previous findings, where all highly abundant miRNAs that show significant changes (Table 1) have been described in previous publications [20,[40], [41], [42], [43], [44], [45], [46], [47], [48], [49]]. The biological implications of the altered miRNAs are consistent with roles in promoting infection. For example, induction of miR-21-5p was found to be driven by the Epstein–Barr nuclear antigen (EBNA)2 [51] and its functional consequence was linked to the deregulation of miR-21-5p target pAKT—an important gene involved in cell proliferation [52]. Induction of miR-146a-5p was linked to the EBV latent membrane protein (LMP)1 via the induction of NFkappaB [47], where it functions in a likely negative feedback loop by modulating the host cell interferon response [42]. In contrast, suppression of miR-200c was previously observed in EBV infected B cells, where it is proposed to maintain latency—as miR-200 family miRNAs promote lytic reactivation [49].

Of the non-micro small ncRNAs, this analysis also recapitulated previous findings for vtRNA1, where the very large (>~600X) fold-change found for vtRNA1 via RT-qPCR (Fig. 2) is consistent with previous results, which found a ~1200-fold increase via subtractive hybridization [50]. It was speculated that the massive induction of vtRNA1 might be to modulate the anti-viral response. This is corroborated by subsequent analyses of other viruses, where vtRNAs were also found to be highly induced: in influenza A virus (IAV), for example, vtRNAs were found to promote viral replication through repressing the activation of PKR and the subsequent antiviral interferon response [53].

In addition to corroborating previous results in additional cell lines, this current study also identified novel non-micro small ncRNAs that are induced by EBV infection. After vtRNA1 the next most consistently highly up-regulated ncRNA was SNORD89, which has not (to our knowledge) been previously identified as being affected by EBV infection. SNORD89, however, was previously found to be up-regulated in IAV-infected cells, alongside many other small ncRNAs [54]. snoRNAs are suggested to play proviral roles in RNA virus infection [55] due to direct interactions between SNORDs and viral transcripts. SNORDs also play roles in promoting DNA virus infection (e.g. CPV and HSV2 [56]). Additionally, EBV encodes its own SNORD analog: v-snoRNA1, which is proposed to be further processed into small bioactive fragments [57]. Our results indicate that induction of SNORD89 (and perhaps others) may also be playing roles in promoting EBV infection.

In addition to SNORDs, tRNAs were also called as being globally induced in BJAB-B1 cells. tRNAs play critical roles throughout viral infection [55] and are essential components of viral gene expression [58]. tRNAs have been previously shown to be induced by herpes viral infection—likely through the global stimulation of pol III transcription, which is essential for EBV ncRNA expression and is posited to contribute to EBV-associated enhancement of cell growth [59].

In this current study we tested six tRNAs called as being highly up-regulated via small RNA-Seq for their induction in EBV-infected cell lines (Fig. 2). Four were validated, with tRNA-Asn (had the second highest calculated fold-change, after vtRNA1; Table A. 2) being the most consistently high across cell types. Although, RT-qPCR results were unable to recapitulate the high up-regulation of tRNA-iMet calculated from RNA-Seq results, we did observe small decreases in this tRNA across EBV-infected cells vs. BJAB (Fig. 2 and Table A. 2). Deregulation of this tRNA is interesting, as a tRNA-iMet precursor was recently shown to be targeted by miR-34a and tRNA-iMet suppression was shown to attenuate cell proliferation while also inducing cell cycle arrest and apoptosis [60]. Thus, there may be relevance to tRNA deregulation on the oncogenic properties of EBV infection. Although miR-34a did not appear as being significantly deregulated between BJAB-B1 and BJAB cells (Table A. 1), other miRNAs may be targeting this tRNA (and others). Additional work is needed to parse out the interactions/effects responsible for ncRNA deregulation associated with EBV infection and disease, as well as the functional consequences of altered host ncRNA levels.

Conclusions

In this study we analyzed differential expression of small ncRNAs between two commonly used laboratory cell lines: BJAB and BJAB-B1. These results help shed light on the importance of considering the roles of non-micro small ncRNAs in EBV infection and disease.

CRediT authorship contribution statement

Lumbini I. Moss: Data curation, Formal analysis, Methodology, Resources, Validation, Visualization, Writing - original draft. Van S. Tompkins: Investigation, Methodology, Writing - review & editing. Walter N. Moss: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing - original draft.