Immune stimuli shape the small non-coding transcriptome of extracellular vesicles released by dendritic cells

The release and uptake of nano-sized extracellular vesicles (EV) is a highly conserved means of intercellular communication. The molecular composition of EV, and thereby their signaling function to target cells, is regulated by cellular activation and differentiation stimuli. EV are regarded as snapshots of cells and are, therefore, in the limelight as biomarkers for disease. Although research on EV-associated RNA has predominantly focused on microRNAs, the transcriptome of EV consists of multiple classes of small non-coding RNAs with potential gene-regulatory functions. It is not known whether environmental cues imposed on cells induce specific changes in a broad range of EV-associated RNA classes. Here, we investigated whether immune-activating or -suppressing stimuli imposed on primary dendritic cells affected the release of various small non-coding RNAs via EV. The small RNA transcriptomes of highly pure EV populations free from ribonucleoprotein particles were analyzed by RNA sequencing and RT-qPCR. Immune stimulus-specific changes were found in the miRNA, snoRNA, and Y-RNA content of EV from dendritic cells, whereas tRNA and snRNA levels were much less affected. Only part of the changes in EV-RNA content reflected changes in cellular RNA, which urges caution in interpreting EV as snapshots of cells. By comprehensive analysis of RNA obtained from highly purified EV, we demonstrate that multiple RNA classes contribute to genetic messages conveyed via EV. The identification of multiple RNA classes that display cell stimulation-dependent association with EV is the prelude to unraveling the function and biomarker potential of these EV-RNAs. Electronic supplementary material The online version of this article (10.1007/s00018-018-2842-8) contains supplementary material, which is available to authorized users.


Introduction
Extracellular vesicles (EV) released by cells are considered as important mediators of intercellular communication [1,2]. These 50-200 nm-sized vesicles are released by a broad range of cells, and have been detected in a wide range of body fluids, such as milk, plasma, urine, and semen [1,3]. The collective term 'EV' refers to a heterogeneous population of secreted vesicles that are formed via different pathways. Exosomes are formed as intraluminal vesicles (ILVs) inside multi-vesicular endosomes (MVE) and are released upon fusion of the MVE with the plasma membrane, whereas microvesicles directly bud off from the plasma membrane [1]. These vesicle populations overlap in size and molecular composition, which currently hampers their discrimination based on biophysical or biochemical parameters [4][5][6].
EV are multi-component entities that transfer proteins, lipids, and RNA between cells. The RNA associated with EV mainly consists of small non-coding RNA species (ncRNA), which are thought to be of major importance in altering molecular processes in the recipient cell [3,4,7,8]. In addition, the RNA composition of EV can give information about the (patho)physiological status of the EV-producing cell. Since EV can be easily obtained from body fluids, such EV-associated RNAs may be used as non-invasive biomarkers for the early detection of disease [9,10]. The molecular composition of EV is often assumed to be a 'snapshot' of the producing cell, from which the tissue type or mutational status of the parental cell can be deduced [9,11,12].
Although the EV-RNA field receives major attention, with quickly rising numbers of publications on this topic, the technical difficulties to obtain pure EV remain underexposed. It is important to note that a substantial amount of extracellular RNA is not contained in EV, but is associated with ribonucleoprotein particles (RNPs), including argonaute 2 (AGO2), or lipoprotein complexes [13][14][15]. These contaminating structures have been shown to coisolate together with EV via commonly used isolation procedures, and may, therefore, greatly affect the outcomes of experiments [4,5,16]. The currently most reliable method to obtain pure EV is by high-speed ultracentrifugation followed by density centrifugation in which EV are separated from contaminating structures based on differences in buoyant density [16]. Subsequently, analyzing the presence of common EV markers (such as CD9, CD63, and CD81) and absence of contaminant proteins should be performed to assess the purity of EV preparations used in the study [17,18].
MiRNAs are the most intensely studied type of small RNA in EV, likely due to the interest in gene-regulatory functions of these RNAs [19][20][21]. It has been shown that changes in the activation status of the EV-producing cell lead to changes in the EV-miRNA composition [22][23][24]. Moreover, miRNAs in EV were shown to be functionally transferred between cells, leading to repression of genes in target cells [25,26]. It is, therefore, thought that differences in the miRNA content of EV underlie distinct functional effects of EV released by differentially stimulated cells [22][23][24]. In addition, changes in EV-miRNA levels in plasma or serum have been associated with diseases such as cancer, rheumatoid arthritis, and Alzheimer's disease [27][28][29].
Despite the interest in miRNAs, a larger part of EVassociated RNA consists of other ncRNA classes. We and others previously showed that EV released by cultured cells and EV present in a variety of body fluids contain many other small non-coding RNA species (20-300 nt), such as tRNA, snoRNA, snRNA, Y-RNA, Vault RNA, and SRP-RNA (also named 7SL) [23,[30][31][32][33][34][35]. Some of these noncoding RNAs are relatively enriched in EV, suggesting that these RNAs are specifically shuttled into EV for release into the extracellular milieu. Inside cells, the above-mentioned non-coding RNAs are known to play a role in basic cellular processes such as translation, RNA splicing, and RNA quality surveillance. Recent studies showed that full-length and fragmented forms of these non-coding RNAs can additionally be involved in other processes including immunological signaling, gene regulation, and guiding of other RNAs and nucleases [36][37][38][39][40][41]. However, it is not known whether activation stimuli imposed on the EV-producing cell affect the levels of these ncRNA classes in EV, similar to EV-associated miRNAs. Evaluation of cell stimulation-dependent changes in EV-associated ncRNA classes is a first step in uncovering their function in EV-mediated signaling processes and their potential as biomarkers indicative of the activation status of cells.
Here, we studied how changes in the activation/differentiation status of the parental cell affect the small non-coding RNA content of EV by employing a primary immune cell model, two different types of cell stimulation, and methods yielding highly pure EV populations. For our analysis, we used primary dendritic cells (DC), which are master regulators of immune responses and have been shown to communicate with various other immune cells via EV [2,8,42]. In addition, the use of DC-derived EV has been proposed as a strategy for cancer immune(chemo)therapy [43][44][45][46]. Depending on external stimuli, DC can be differentiated into functionally different phenotypes that either activate or down-regulate immune responses [47]. These DC, therefore, present a suitable model to investigate how different external stimuli affect incorporation of a broad range of RNA classes into EV. We compared cellular and EV-associated RNA levels in control, immune-stimulating, or immune-suppressing DC conditions using RNA deep sequencing and RT-qPCR. Our data indicate that differentiation signals imposed on dendritic cells affect the EV-associated levels of particular RNA classes, such as miRNA, snoRNA, and Y-RNA, while other non-coding RNA types remain largely unchanged. Moreover, only a minor part of the stimulation-induced changes in EV-RNA content reflects changes in cellular RNA levels. This study exemplifies how comprehensive analysis of RNA obtained from highly purified EV yields candidate small ncRNAs beyond miRNAs for further exploration as biomarkers or functional entities within EV.

Cell culture
Complete culture medium was prepared by supplementing Iscove's Modified Dulbecco's Medium (IMDM, Lonza, Verviers, Belgium) with 2 mM Ultraglutamine (Lonza, Verviers, Belgium), 10% fetal calf serum (FCS, GE Healthcare Bio-Sciences, Pasching, Austria), 100 IU/ml penicillin and 100 μg/ml streptomycin (Gibco, Paisley, United Kingdom), and 50 µM β-mercaptoethanol. The GM-CSF producing NIH 3T3 cell line (R1) was grown in complete medium. Primary bone marrow cells were grown in complete medium supplemented with 30% conditioned medium from R1 cells. To prepare EV-depleted medium, a mixture of 150 ml conditioned R1 culture supernatant and 42.5 ml FCS (end concentration 30% FCS) was depleted of EV by overnight centrifugation at 100,000g in an SW28 rotor (k-factor 334.2) (Beckman Coulter, Brea, CA). The EV-depleted supernatant was carefully pipetted from the tubes, leaving 5 ml in the tubes to prevent disturbance of the pellet. This supernatant was filtered through a 0.22 μm bottle top filter (Millipore, Billerica, MA) after which additional IMDM, ultraglutamine, antibiotics, and β-mercaptoethanol were added to prepare complete medium as indicated above. The efficiency of EV depletion from culture medium, as regularly assessed by high-resolution flow cytometric analysis, is on average ~ 90%. Primary bone marrow cells were flushed from the femur and tibia of 8-12-week-old C57bl/6 mice and differentiated to dendritic cells according to Lutz et al. [48]. To generate tolerogenic DC, cells were treated from day 2 onward with 10 nM 1α,25-dihydroxyvitaminD3 (Sigma, St Louis, MO). To generate immunogenic DC, 1 µg/ ml lipopolysaccharide (LPS, cat. L-2630, Sigma, St. Louis, MO) was added on day 12. On day 13, non-adherent and semi-adherent cells were recovered and cultured at 3 × 10 6 cells/dish for 20 h in EV-depleted culture medium. Cell viability, as determined by Trypan blue exclusion, did not differ between treatments and was above 90% for all cultures. All cells were maintained at 37 °C and 5% CO 2 in a humidified incubator. Experiments were approved by the institutional ethical animal committee at Utrecht University (Utrecht, The Netherlands).

Fluorescent labeling and purification of EV
Conditioned cell supernatants were pooled to volumes of 100-140 ml per condition and were subjected to differential centrifugation as described previously [49]. In brief: supernatant was sequentially centrifuged 2 × 200g for 10 min, 2 × 500g for 10 min, and 1 × 10,000g for 30 min. Next, EV were pelleted by ultracentrifugation at 100,000g for 65 min using an SW28 rotor (k-factor 334.2) (Beckman Coulter, Brea, CA). For EV quantification by high-resolution flow cytometry, pellets were resuspended in 20 µl PBS + 0.2% BSA (cleared from aggregates by overnight ultracentrifugation at 100,000g) and labeled with 1.5 µl PKH67 (Sigma, St. Louis, MO) in 100 µl Diluent C per pellet. For EV-RNA isolation, 100,000 g pellets were resuspended in 50 µl PBS + 0.2% BSA. PKH67-labeled samples or samples for RNA isolation were mixed with 1.5 ml 2.5 M sucrose, and overlaid with a linear sucrose gradient (2.0-0.4 M sucrose in PBS). Gradients were centrifuged 15-18 h at 192,000g in a SW40 rotor (k-factor 144.5) (Beckman Coulter, Brea, CA). PKH67-labeled EV were used for high-resolution flow cytometric analysis (see below). Alternatively, the fractions with densities of 1.12-1.18 g/ml (as measured by refractometry) were pooled, diluted six times in PBS + 0.2% EV-depleted BSA, and centrifuged again at 192,000g for 65 min in a SW40 rotor (k-factor 144.5) prior to EV-RNA isolation.

RNA isolation
Small RNA was isolated from EV pellets and from 1 × 10 6 cells using the miRNeasy Micro kit according to the small RNA enrichment protocol provided by the manufacturer (Qiagen, Hilden, Germany). RNA yield and size profile were assessed using the Agilent 2100 Bioanalyzer with Pico 6000 RNA chips (Agilent Technologies, Waldbronn, Germany).

Preparation of small RNA sequencing libraries
50 ng cellular small RNA and 2 ng EV-derived small RNA was treated with DNase [Turbo DNA-free kit (Life Technologies, Carlsbad, CA)] according to the manufacturer's instructions. The DNase-treated RNA was pelleted using Pellet Paint (Merck, Darmstadt, Germany) according to the manufacturer's instructions, reconstituted in milliQ (MQ) and subsequently subjected to ribosomal RNA depletion using RiboZero Gold kit (Human/Mouse/Rat) (Illumina, San Diego, CA) according to the manufacturer's instructions. This is required to deplete residual rRNA from cellular RNA samples and, for comparability, we subjected EV-RNA samples to the same procedure. Hereafter, RNA was pelleted with Pellet Paint and reconstituted into 6 µl MQ. cDNA libraries were prepared using the NebNext smallRNA library prep kit for Illumina (New England Biolabs, Ipswich, MA), according to the manufacturer's instructions but with the following adaptations: 3′ adapter-ligation was carried out overnight at 16 °C; PCR amplification was done using Kapa HiFi Readymix 2× PCR mastermix (Kapa Biosystems, Wilmington, MA) using barcoded primers and the following PCR programme: 2 min at 95 °C, 15 cycles of 20 s at 98 °C, 30 s at 62 °C, 15 s at 70 °C, and a final elongation step of 5 min at 70 °C. cDNA was purified using magnetic AMPure XP beads (Beckman Coulter, Brea, CA) and quantified using Agilent 2100 Bioanalyzer and DNA HiSensitivity chips (Agilent Technologies, Waldbronn, Germany). Adapter dimers (126 nt in size) were removed by running the cDNA on 6% TBE gels (Life Technologies, Carlsbad, CA) for 60 min at 145 V, after which cDNA products of 15-300 nt (+ 126 nt adapters) were cut from the gel and purified. Subsequently, all libraries were pooled at equimolar ratios and run on a 4-12% TBE gel (Life Technologies, Carlsbad, CA). Size fractions of 15-25 nt (includes miR-NAs), 25-60 nt, 60-80 nt (includes tRNAs), and 80-275 nt (each + 126 nt adapters) were cut from the gel, and purified and pooled in a volumetric ratio of 2:2:3:3. Sequencing was done using 150 bp paired-end reads on an Illumina HiSeq 4000 machine (Illumina, San Diego, CA) at ServiceXS (Leiden, The Netherlands).

RNAseq data analysis
All data was processed as paired-end reads. Data quality was checked with FastQC and reads were processed with cutadapt (version 1.8) to remove low-quality reads, clip the reads to 100 bp, and remove adapter sequences.
Differential abundance was determined separately for each ncRNA class, and separately for EV and cells, using the edgeR package in R [50]. Data were normalized using the TMM method (weighted trimmed mean of M values). Estimating the common, trended (overexpression values), and tagwise dispersion, a generalized linear model, containing the effect of cell stimulation and the effect of the day-to-day variation between triplicate experiments, was fit. The log-likelihood ratio test was used to evaluate differential expression between treatments relative to control. p values were adjusted for multiple testing using Benjamini and Hochberg's false discovery rate (FDR). Average fold-change over three independent experiments and standard deviation were plotted. Analysis of RNA fragments was done using the UCSC genome browser and Integrated Genome Viewer [51].

Quantitative real-time PCR
cDNA was generated from cellular or EV-derived small RNA using the miScript RT2 kit (Qiagen, Hilden, Germany). An equivalent of 20 pg RNA was used per qPCR reaction and mixed with 100 nM primers (Isogen Life Sciences, De Meern, The Netherlands) and 4 µl SYBR Green Sensimix (Bioline Reagents Ltd., United Kingdom) in an 8 µl reaction. No-RT-controls confirmed the absence of genomic DNA and non-specific amplification.
Cycling conditions were 95 °C for 10 min followed by 50 cycles of 95 °C for 10 s, 57 °C for 30 s, and 72 °C for 20 s. All PCR reactions were performed on the Bio-Rad iQ5 Multicolor Real-Time PCR Detection System (Bio-Rad, Hercules, CA). Quantification cycle (Cq) values were determined using Bio-Rad CFX software using automatic baseline settings. Thresholds were set in the linear phase of the amplification curve.

High-resolution flow cytometric analysis of EV
High-resolution flow cytometric analysis of PKH67labeled EV was performed using a BD Influx flow cytometer (BD Biosciences, San Jose, CA) with an optimized configuration, as previously described [49,52]. In brief, we applied threshold triggering on fluorescence derived from PKH67-labeled EV passing the first laser. Forward scatter (FSC) was detected with a collection angle of 15°-25° (reduced wide-angle FSC). Fluorescent 100and 200-nm polystyrene beads (FluoSpheres, Invitrogen, Carlsbad, CA) were used to calibrate the fluorescence and rw-FSC settings. Sucrose gradient fractions containing PKH67-labeled EV were diluted 25× in PBS and vortexed just before measurement. This dilution factor was sufficient to avoid 'coincidence' (multiple EV arriving at the measuring spot at the same time), thereby allowing accurate quantitative comparison of EV numbers in different conditions. Moreover, samples were measured at maximally 10,000 events per second, which is far below the limit in the electronic pulse processing speed of the BD Influx [53].

Nanoparticle tracking analysis (NTA)
EV purified by density gradient ultracentrifugation as described above were resuspended in 50 µl PBS + 0.4% EV-depleted BSA. Samples were diluted 20 times in PBS (confirmed to be particle-free when analyzed with the same settings as used for the EV samples) before measurement on a Nanosight NS500 instrument (Malvern, Worchestershire, UK). Data acquisition and processing were performed using NTA software 3.1. Each sample was recorded for 5 times 30 s, 25 frames per second, camera level 14, detection threshold 5.

Electron microscopy
EV and RNP were separated by density gradient centrifugation as described above, and resuspended in 15 µl PBS. 3 µl aliquots of either EV or RNP were added on top of carboncoated 300 mesh, copper grids, and incubated for 2 min at room temperature. The grids were then washed two times in 50 µl PBS and stained with 2% uranyl acetate. Imaging was performed on a T20 electron microscope (FEI) operated at 200 keV. Images were recorded on a CCD Eagle camera (FEI).

Northern blotting
Cellular or EV-derived small RNA (between 100 ng and 1 ng per lane, as indicated) was denatured in gel loading buffer II (Invitrogen) for 2 min at 70 °C and snap-cooled on ice before loading on a denaturing 15% PAGE gel (National Diagnostics, Nottingham, UK) as described previously [54]. RNA was transferred onto a Hybond-N Nylon membrane (Amersham Pharmacia Biotech, GE Healthcare Bio-Sciences, Pasching, Austria) and chemically crosslinked with 1-ethyl-3-(3-dimethylaminopropyl) carbodimide (Sigma) at 55 °C for 2 h [55]. 32 P-labeled DNA oligo probes perfectly complementary to each RNA species were hybridized overnight at 42 °C in PerfectHyb (Sigma) solution. Blots were analyzed by phosphorimaging using a Typhoon Scanner (GE Healthcare).

RT-qPCR probes
The following forward primers were used in RT-qPCR (5′ to 3′) in combination with the miScript universal reverse primer.

Statistical analyses
Differences between the numbers of EV released by LPSand VitD3-stimulated DC versus control DC were analyzed by one-way ANOVA with Dunnett's two-sided post hoc test. Similar statistical testing was performed on differences in RNA yield. Differences in RNA biotypes in EV and cells were analyzed by two-tailed t test. Differences in foldchanges of individual RNAs as measured in RT-qPCR were analyzed by one-way ANOVA. Significance was defined as p < 0.05. Statistical tests on data not derived from RNA sequencing were done in SPSS (v24, IBM).

Availability of data and materials
Raw data and count tables are deposited in the GEO database under accession number GSE105151. Written details on experimental procedures have been submitted to the EV-TRACK knowledgebase (EV-TRACK ID: EV170030) [5].

Results
Functionally distinct subsets of mouse bone marrowderived DCs were generated by exposure of cells to LPS, thereby inducing immune-stimulatory DC (LPS-DC), or to 1α,25-dihydroxyvitaminD3, thereby inducing immunesuppressive DC (VitD3-DC) [56]. Untreated DC were used as control cells. The different DC subtypes displayed the expected differences in expression of the activation markers, with LPS-DC showing increased levels of MHC class II and CD86 expression, while VitD3-DC showed increased PD-L1/CD86 ratios [57] and resistance to LPS activation [58] (Supplementary Fig. 1A-F). We isolated EV from cell culture supernatant of LPS-DC, VitD3-DC, and control DC using differential (ultra)centrifugation followed by density gradient separation, as previously published by our group (Supplementary Fig. 2A and [52]). Transmission electron microscopic analysis indicated that low-density fractions contained 100-150 nm-sized EV, whereas no EV were observed in high-density fractions enriched in protein complexes ( Supplementary Fig. 2B). Our in-house developed high-resolution flow cytometric method was used for * * n.s. high-throughput EV analysis at the single particle level [49,52]. The light scatter patterns induced by EV from differently stimulated DC were highly similar ( Supplementary  Fig. 2C). However, we detected differences in the number of released EV, with LPS-DC releasing more, and VitD3-DC releasing less EV than non-treated control DC (Fig. 1a). These quantitative differences were confirmed by nanoparticle tracking analysis (NTA) (Supplementary Fig. 2D). In addition, the NTA data indicated that the majority of EV observed in all conditions were in the 100-200 nm size range. Western blot analysis showed that the DC-derived EV contained variable levels of common EV proteins such as CD9, MHCII, CD63, and Galectin-3 [59], whereas abundant cellular proteins such as beta-actin were not detected ( Supplementary Fig. 2E). These data demonstrate that VitD3 and LPS treatment of DC differentially affect the number and protein content of EV released by these cells. Next, we analyzed how LPS-or VitD3 treatment of DC affected the RNA content of EV released by these cells. It is important to note that extracellular RNA can be associated with either EV or ribonucleoprotein complexes, which may sediment at similar centrifugal force [4]. Using density gradient-based purification, we separated EV from contaminating ribonucleoprotein complexes. Depending on the differentiation status of the parental DC, the EV-free ribonucleoprotein fraction contained between 26 and 55% of total extracellular RNA released from cells (Supplementary Fig. 2F-G). The total amount of EV-associated RNA released by equal numbers of control, LPS-, and VitD3treated DC was different and reflected the differences in EV numbers (Fig. 1a, b). To assess whether DC stimulation altered the RNA composition of EV, we prepared sequencing libraries of EV-associated and cellular small RNA from control-, LPS-, and VitD3-DC cultures (n = 3 biological replicates). Sequencing libraries were prepared using an adapter-ligation-based method routinely used in miRNA profiling [60]. Such a method has been frequently applied in studies to demonstrate the presence of miRNA and other small non-coding RNA types in EV (e.g., [23,24,31,32,61]). All libraries contained RNAs with lengths ranging between 20 and 300 nucleotides. Since short RNA molecules have a selective advantage during PCR amplification and pre-amplification on the sequencing flow-cell, leading to an apparent over-representation of short read lengths, we enriched the libraries for longer RNAs to ensure sufficient coverage of these sequences. To this end, cellular and EV-associated cDNA libraries were split into different size fractions of 15-25 nt (includes miRNAs), 25-60 nt, 60-80 nt (includes tRNAs), and 80-275 nt. Subsequently, we enriched for longer RNA molecules by pooling these size fractions in a volumetric ratio of 2:2:3:3 (for details, see "Materials and methods"). This, indeed, resulted in a better coverage of 100-300 nt mid-size RNAs, such as full-length Vault RNA (142 nt) (Supplementary Fig. 3). More than 13 million reads were obtained for each library (Supplementary Table 1). Pearson coefficients for biological replicate samples were ≥ 0.85 for EV, and ≥ 0.95 for cells (Supplementary Fig. 4). As expected based on the previous studies [23,[30][31][32]34], we detected the presence of various small RNA classes in EV (miRNA, snoRNA, snRNA, tRNA, mis-cRNA, and rRNA). Some RNA types were enriched in EV, such as Y-RNA, Vault RNA, and 7SL RNA classified as 'miscRNA', while other RNA types, such as snoRNA, were relatively less abundant in EV compared to cells (Fig. 2). Since we depleted ribosomal RNA before preparation of the sequencing libraries, apparent enrichment of rRNA in EV is probably caused by differences in depletion efficiency between cells and EV. MiRNA and snRNA abundance was similar between cells and EV. Stimulation of DC with LPS or VitD3 did not lead to major changes in the distribution of sequencing reads over different RNA classes in EV and cells ( Supplementary Fig. 5A-B).
Next, we assessed whether, within each of the different classes, the levels of individual RNAs in EV changed as a result of cell stimulation. Hereto, read counts were normalized per RNA class, after which the fold changes in reads The scatter plots of these data indicated that LPSand VitD3-treatment affected the levels of miRNA, snoRNA, snRNA, and miscRNA in EV (Fig. 3a-e). Some of the stimulation-induced changes were similar for LPS and VitD3 conditions, while other RNAs show opposing alterations, or change only in one of the stimulation conditions. Of the tested RNA classes, tRNA levels in EV showed the lowest rate of change due to LPS or VitD3 stimulation. 97% of tRNA reads mapped to six different tRNA isoacceptors (GluCTC, GlyGCC, LysCTT, GluTTC, GluCTG, and LysTTT) (Supplementary Fig. 6A). Although relatively high tRNA read counts were observed in the EV-associated RNA pool ( Supplementary Fig. 5B), the levels of these six abundant tRNAs remained stable (log2FC < 1) in response to cellular stimulation (Supplementary Fig. 6B). Besides tRNAs, the highly purified EV populations used in this study also contained substantial levels of snRNA, which is an RNA type known to be mainly confined to the nucleus. Although cell stimulation seemed to induce changes in EV-associated snRNA levels (Fig. 3c), many of the individual data points mapped to multicopy genes with highly similar sequences (> 95% sequence identity) corresponding to known snRNAs. Cumulative analysis of these multicopy genes indicated that 96% of snRNA reads mapped to four different snRNAs (U1, U2, U5, and U6) (Fig. 4a), of which the levels in EV did not significantly differ between stimulation versus control conditions (Fig. 4b). Thus, the levels of snRNAs and tRNAs, Fig. 3 LPS-versus VitD3-induced changes in EV-associated RNA classes. EV-associated RNA from control, LPS, and VitD3 conditions was isolated and analyzed by RNA sequencing. Read counts for individual RNAs were normalized to the total read counts of each RNA class. a-e LPS-or VitD3-induced fold changes and corresponding p values were calculated relative to the control condition with edgeR GLM method. Data are expressed as log2fold change in LPS-EV relative to control-EV versus log2fold change in VitD3-EV relative to control-EV. Cutoffs for log2fold changes larger or smaller than 1 are indicated with dashed lines (red = LPS; green = VitD3), so all data points beyond these lines are differentially expressed with log2FC > 1. Dot size represents the normalized abundance (logCPM) of individual RNAs. RNAs that changed with non-adjusted p < 0.05 are highlighted in black which are commonly detected RNA types in EV, remained constant upon different immune stimuli imposed on the EVproducing cells.
In contrast to tRNAs and snRNAs, the EV-associated levels of three other RNA types tested here, i.e., miRNAs, snoRNAs, and miscRNAs, changed in response to the different stimulate imposed on the DC. We aimed to identify the RNAs differentially present in EV released by LPS-versus VitD3-treated DC, since such RNAs may underlie distinct effects of EV on recipient cells.
MiRNAs known to be involved in the pro-inflammatory function of DC were enriched in LPS-induced EV, whereas EV from VitD3-treated DC were enriched in miRNAs that dampen immune-stimulatory signaling cascades in DC. We validated three LPS-enriched and three VitD3-enriched RNAs by RT-qPCR. The use of stable reference transcripts in RT-qPCR analysis is important to ensure reliable normalization between different samples. However, the identity and general applicability of reference transcripts stably present in EV are understudied. Here, we used our sequencing data set to select four non-coding RNAs that were present in EV at comparable levels across all conditions, and confirmed their stability by RT-qPCR ( Supplementary Fig. 7A-B). We normalized our qPCR data to the geometric mean of these four genes to minimize errors caused by technical and/or biological variation [74]. Using this normalization strategy, RT-qPCR analysis confirmed the sequencing data for the tested miRNAs, although miR146-5p did not reach significance (Fig. 5b). Together, these data indicate that opposing stimuli imposed on the same type of parental cell induce the release of functionally different sets of miRNAs via EV.
Besides changes in miRNA levels, the levels of several EV-associated snoRNAs were found to differ between LPS and VitD3 conditions (Fig. 5c), with three H/ACA box snoRNAs (snora2b, snora32, and snora55) displaying FDR levels < 0.05. The overall low abundance of this RNA class in EV (Fig. 2c) hampered reliable quantification of snoR-NAs in EV by RT-qPCR (data not shown). Nevertheless, our data are a first indication that EV-associated snoRNA levels change with the activation status of parental cells. Although none of the RNAs in the miscRNA group reached FDR levels < 0.05, analysis by RT-qPCR indicated that the levels of Y3-RNA (Rny3) and an additional member of the Y-RNA family, Y1-RNA (Rny1), significantly differed between LPS and VitD3 conditions (Fig. 5d, e). In cells, the conserved family of Y-RNAs is known to bind to several different proteins that play a role in RNA quality control [75,76]. Members of the small non-coding Y-RNA family have frequently been detected in EV from multiple different cell types and different body fluids [30-32, 34, 35, 77-81]. We here provide evidence that Y-RNA levels in EV can be regulated by stimuli imposed on the EV-producing cell. The reduced Y-RNA levels in EV from LPS-stimulated DC compared to EV from VitD3-stimulated DC (Fig. 5e) indicate that the presence of Y-RNAs in EV may be indicative for the immune function of the parental cells. Overall, our data indicate that immune cell stimuli cause changes in the EV-associated levels of specific miRNAs, snoRNAs and Y-RNAs.
So far, we compared the abundance of specific RNA classes in EV by assessing the total numbers of reads mapping to non-coding RNA genes. However, also the  32,34,35,77,78,81]. Especially, Y-RNA fragments (19-35 nt) have gained attention, because they were found to associate with Argonaute, suggesting scope for gene-regulatory functions [82,83], and because their presence in plasma has been linked to disease [78,84]. We examined our data to determine whether reads mapping to Y-RNAs included these fragments in the purified DCderived EV. Based on RNAseq data, 5′ and 3′ fragments of Y1-RNA were present in EV and such fragments seemed more abundant than the full-length form of this Y-RNA (Fig. 6a). To confirm the full-length or fragmented nature of the Y-RNA molecules, we performed Northern blot analysis. Using probes recognizing the 5′ and 3′ end of Y1-RNA, we found that both 5′ -and 3′ fragments and full-length Y-RNA could be detected in both cellular and EV-RNA (Fig. 6b-d).
Strikingly, the relative amount of full-length Y1-RNA detected by Northern blot was much higher than expected based on the sequencing coverage plots. This indicates that RNAseq-based detection of the full-length form of Y1-RNA was highly inefficient. This is likely due to the 5′-triphosphate group of full-length Y-RNAs or their strong secondary structure, which may hamper sequencing adaptor ligation and thereby efficient amplification and sequencing of these RNAs. Interestingly, Northern blot analysis also indicated that the levels of full-length Y1-RNA in LPS-EV were reduced, which was in accordance to the RT-qPCR analysis for full-length Y1 (Fig. 5e), while the levels of fragmented Y1-RNA in EV remained relatively stable (Fig. 6d). These data urge caution in classifying fragmented and full-length forms of Y-RNA based on RNA sequencing data. Changes in the RNA content of EV have been suggested to directly reflect changes in the RNA levels of the parental cell. To test this hypothesis, we investigated whether all of the stimulation-induced changes in miRNAs, snoR-NAs, and miscRNAs observed in EV reflected changes in the RNA pool of the parental cells and vice versa. After normalization within each RNA class, fold-changes in EV were plotted against the fold changes in cells for miRNA, miscRNA, and snoRNA (Fig. 7a), with dashed lines indicating log2FC of 1 or − 1. Dot size indicates the average (n = 3) abundance of a transcript in EV, and dot coloring is used to indicate RNAs that were significantly changed in cells, EV, or both. First, we tested whether changes in EV-associated RNAs reflected stimulationinduced changes in cells. RNAs of which the levels in EV significantly changed (non-adjusted p values < 0.05) upon LPS or VitD3-treatment were categorized as 'reflective RNAs' (p < 0.05 in both cells and EV), 'EV-responsive RNAs' (p < 0.05 in EV, p > 0.05 in cells), or 'reciprocal RNAs' (opposite changes in EV and cells with p < 0.05) (Supplementary Table 2). Of all RNAs that significantly changed in EV (set to 100%), a minor percentage of the miRNAs and snoRNAs reflected the up/downregulation in cells or were reciprocally regulated, while most of the RNAs changed only in EV, but not (significantly) in cells (Fig. 7b-d). Conversely, we aimed to assess which proportion of LPS-or VitD3-induced changes in cellular RNA was reflected in the EV-RNA content. RNAs that changed in cells upon cell stimulation (non-adjusted p values < 0.05) were, therefore, classified as 'reflective' (p < 0.05 in both cells and EV), 'cell-responsive' (p < 0.05 in cells, p > 0.05 in EV), 'reciprocal' (opposite changes in EV and cells with p < 0.05), or 'not in EV' (p < 0.05 in cells, but not detected in EV). Comparable to what we observed in the above analysis of EV-RNA, only a minor percentage of changes in cellular miRNAs and snoRNAs was reflected in EV. The majority of the RNAs that changed in cells did not significantly change in EV (Fig. 7e-g), and 10-20% of the RNAs that changed in cells were not detected in EV. We selected three reflective and three EV-responsive RNAs (Fig. 7h) for validation by RT-qPCR. For the reflective miRNAs miR-155, -9, and -10a, we confirmed that the stimulation-induced fold-change in EV-associated and cellular RNA was highly similar. For the EV-responsive miRNA miR-146a and the two Y-RNAs, significant cell stimulation-induced changes were observed in EV, while cellular levels did not change significantly compared to the control conditions. Importantly, these data indicate that part of the cell stimulationinduced changes in EV-RNA content match the changes observed in cellular RNA, but that cell stimulation also induces changes in RNA levels that are only observed in cells or in EV.
Overall, the results obtained for our primary DC cultures show that the molecular messages enclosed in EV are composed of multiple RNA classes that are specifically shuttled into EV dependent on the activation status of the Fig. 5 EV from LPS-/VitD3-stimulated DC display differences in the levels of several miRNAs, snoRNAs, and Y-RNAs. EV-associated RNA from control, LPS, and VitD3 conditions was isolated and analyzed by RNA sequencing. Read counts for individual RNAs were normalized to the total read counts of each RNA class. LPS-or VitD3-induced fold changes and corresponding p values were calculated relative to the control condition. a Volcano plots of EV-associated miRNAs in LPS versus VitD3 conditions. Thresholds for twofold change and non-adjusted p < 0.05 are indicated. Data points represent average values of n = 3 biological replicates. Significant changes are indicated with different colours. Grey: non-significant, blue: nonadjusted p value < 0.05, and red: FDR < 0.05. b RT-qPCR validation of six miRNAs showing differential abundance between LPS-EV and VitD3-EV. Data are expressed as log2fold change in LPS-and VitD3-EV compared to control-EV. Indicated are the mean ± SD values of n = 4 independent experiments, one-way ANOVA, *p < 0.05. Volcano plots of EV-associated snoRNAs (c) and miscRNAs (d) in LPS versus VitD3 conditions. Thresholds for twofold change and non-adjusted p < 0.05 are indicated. Data points represent average values of n = 3 biological replicates. Significant changes are indicated in different colours. Grey: non-significant, blue: non-adjusted p value < 0.05, red: FDR < 0.05. e RT-qPCR validation of Y3 and its family member Y1 showing differential abundance between LPS-EV and VitD3-EV compared to control-EV. Indicated are the mean ± SD values of n = 4 independent experiments, one-way ANOVA, *p < 0.05 ◂ Table 1 MicroRNAs enriched in LPS-or VitD3-EV with known functions in DC EV-associated RNA from control, LPS, and VitD3 conditions was isolated and analyzed by RNA sequencing. Read counts for individual RNAs were normalized to the total read counts of each RNA class. LPS-or VitD3-induced fold-changes and corresponding p values were calculated relative to the control condition with edgeR GLM method. We created a top 20

Discussion
The data presented here broaden our view on the plasticity of the small RNA content of EV in relation to changes in the activation status of the EV-producing cell. Whereas most of the previous studies focused on the miRNA content of EV, we show that cell status-dependent changes in the RNA composition of EV also extend to other small RNA classes, which may, therefore, also contribute to the specific genetic messages conveyed by EV to recipient cells. The immune cell model which we employed, i.e., DC that are well known for their strong and diverse functional responsiveness to external stimuli, served two purposes. First, the cell system illustrates that diverse stimuli can induce the release of EV with variable levels of multiple RNA types, and that only some of the stimulation-induced changes in EV-associated RNAs reflect changes in the cellular RNA. Second, the data extend our knowledge on the central role of DC in raising and regulating immune responses. The identification of multiple EV-RNA types that are specifically associated with either the immuneactivating or immune-suppressing status of DC is the prelude to unraveling the function of these RNAs in EV.
So far, RNA sequencing studies in which multiple RNA classes were identified in EV employed EV isolation methods known to co-isolate extracellular protein-RNA complexes (RNP) that overlap in size with EV, e.g., ultracentrifugation and precipitation-based methods [4,24,30,32,34,80,85,86]. The results of our present study indicate that up to 55% of total extracellular RNA present in ultracentrifugation pellets of cell culture supernatant can be associated with RNPs. Buoyant density-based separation of EV from RNP, as employed in this study, allowed us to specifically identify RNA species that cells differentially sort into EV upon exogenous stimulation. In RNA sequencing libraries prepared from these highly pure EV, we observed that the percentage of alignment to the mouse genome was lower than in cellular RNA libraries. This is likely caused by the (ultra)low-input quantity of RNA in our EV-sequencing libraries, which amplifies the contribution of lab-derived contaminant RNAs, e.g., from commercial nucleic acid extraction kits or sample cross-contamination [87,88]. The percentage of alignment may also be dependent on the size range of RNAs selected for sequencing. Whereas, in most published EV-RNA sequencing studies, RNAs in the 20-40 nt size range were sequenced, thereby strongly selecting for microR-NAs, we here sequenced RNAs of a much larger size range (20-300 nt). We previously showed that several of these larger and highly modified/structured RNAs, for which the sequencing efficiency may be compromised, are enriched in EV compared to cellular RNA [30], which may have caused the percentages of alignment to be lower in EV compared to cells. For our studies, we selected a commercial adapter-ligation-based method that has been frequently used in EV-RNA sequencing studies [23,24,31,32,61]. While optimized for miRNA analysis, it should be noted that these library preparation methods show bias in the efficiency with which sequences with base modifications, terminal triphosphate groups, and strong secondary structures are ligated and sequenced [89]. As an example, we here provide Northern blot-based evidence that the frequently reported predominance of Y-RNA fragments over full-length Y-RNA in EV [30,32,35,90] represents a sequencing artefact. We observed a similar discrepancy between sequencing read coverage and Northern blot analysis for tRNAs ( Supplementary Fig. 8). While the triphosphate group on the 5′-end [91] or their strong secondary structure may hamper efficient sequencing of full-length Y-RNA, difficulties with sequencing of full-length tRNAs are likely due to the inability of the reverse transcriptase enzyme used during sequencing library preparation to read through the highly modified tRNA structure [92]. This is corroborated by recent sequencing studies deploying reverse transcriptases that are insensitive to secondary RNA structures [93,94]. Together, these findings urge caution in the interpretation of RNA fragmentation based on RNA sequencing data alone. Moreover, we reduced the effect of sequencing biases on the assessment of quantitative differences in EV-RNA content by determining the fold changes in identical transcripts between different conditions.
Our experimental approach allowed comprehensive analysis of changes in EV-RNA classes induced by diverse immune-relevant stimuli. The most prominent stimulationinduced changes in EV-associated RNA levels were found for miRNA, Y-RNA, and snoRNAs (Fig. 3). Although detected in EV from multiple sources [23,31,[33][34][35], the presence of snoRNAs in EV has been largely understudied. Here, we provide evidence that the levels of EV-associated snoRNAs can be regulated by stimuli imposed on the EV-producing cell. Since changes in the cellular levels of snoRNA have been associated with diseases [36,95,96], our present data point to EV-associated snoRNAs as potential functional entities in EV and interesting candidate biomarkers indicative of the cellular activation status.
Numerous studies already demonstrated that miRNA levels in EV change upon disease induction or cell stimulation [22][23][24]. Unique aspects of our present RNAseq study are the use of two different stimuli for comparative analysis of RNA levels in highly purified EV and the parallel assessment of cellular and EV-associated RNA. The c d e h f g enrichment of endotoxin-responsive miRNAs, such as miR-155 and miR-146a [20,97,98], that we observed in EV from LPS-stimulated DC concurs with previously published semi-quantitative microarray-based data on the miRNA content of EV from stimulated DC [22,26]. Both miR-155 and miR-146a are known to be endotoxin-responsive, but miR-155 has an immune-activating role, while miR-146a is involved in dampening of immune responses and could act as a molecular brake on inflammation [99]. It is thought that the coordinated action of these two is important in regulation of immune response [21]. Differences in the miR-155/ miR-146a ratios that we observed in LPS-EV and VitD3 EV (Fig. 5b) may, therefore, lead to different effects of these EV on immune activation. In addition, we here demonstrated that a different (i.e., immune downregulatory) type of stimulus imposed on DC led to both a reduction of EV-associated levels of these immune-activating miRNAs and an enrichment in miRNAs implicated in the suppression of immune responses via modulation of multiple pathways. Examples include dampening of TLR signaling by miR-27a and miR-126a [70,73], and downregulation of (regulators of) PI3K/ Akt and NF-κB by miR-378 [100], miR-708 [71] and miR-27b-3p [101], processes that have been implicated in the suppression of pro-inflammatory cytokine release. Such differences in the EV-miRNA content may therefore contribute to immune-activating or immune-suppressive functions of differentially stimulated DC. Not only immune-suppressive miRNAs but also two types of Y-RNA were found to be specifically excluded from EV released by LPS-stimulated DC (Fig. 5e). This raises the question whether Y-RNAs contribute to immune-related functions of EV. Until now, different types of Y-RNA (mY1 and mY3 in mice and hY1, hY3, hY4, and hY5 in humans) have been detected in EV from a multitude of different cell types and in different body fluids [30-32, 34, 35, 77-81]. The extracellular presence of Y4-and Y5-RNA and fragments thereof has mostly been associated to cancer [32,81,84,102] and coronary artery disease [78,103]. Immunerelated functions of Y-RNAs include their capacity to trigger TLR signaling. Interestingly, Y-RNAs may vary in their specificity for different TLRs, with Y3-RNA triggering predominantly TLR3, while Y1-, Y3-, and Y4-RNA trigger TLR7 [104]. On the contrary, the reduction of Y1-and Y3-RNA levels in EV released by immune-stimulatory DC, shown in this study, suggests an immune downregulatory role for these EV-associated Y-RNAs. In line with this, high levels of Y-RNAs have been reported in immune-suppressive EV released by the parasite Heligmosomoides polygyrus and in seminal fluid EV for which immunosuppressive effects have been described [35,[105][106][107]. Overall, the observed variable presence of Y-RNA in EV from immune cells has raised further interest in unraveling the immune-related function and biomarker potential of these RNAs in EV.
With regard to the biomarker potential of EV-associated RNAs, our data also urge caution in interpreting EV as snapshots of the cell from which they arise. We found that only a subset of changes in EV-RNA reflected changes that occurred at the cellular level. The highest proportion of reflective RNAs was observed within the category of miR-NAs and included prominent immune-related RNAs such as miR-155, miR-9, and miR-10-5p (Fig. 7h). For a substantial number of EV-associated miRNAs and both Y-RNA types, however, stimulation-induced changes were observed in EV, but not in cells. This could suggest that EV-levels of reflective RNAs are mainly regulated by transcription, whereas levels of EV-responsive RNAs are mainly determined by shuttling rate. For RNAs known to primarily reside in the nucleus (such as snRNAs and snoRNAs), subtle changes in the cytoplasmic pool of these RNAs may be overshadowed by the much larger nuclear pool. Comparison of EVassociated levels with cytoplasmic instead of total cellular RNA levels may, therefore, provide a better insight in shuttling of these nuclear RNAs into EV. Interestingly, we also observed that many of the stimulation-induced changes in cellular RNA levels were not reflected in the RNA levels of EV released by these cells. This strengthens the hypothesis that cells release EV with selective sets of RNAs into the extracellular milieu. From a biomarker perspective, these Fig. 7 LPS-/VitD3-induced changes in cellular versus EV-associated RNA levels. Cellular and EV-associated RNA from control, LPS, and VitD3 conditions was isolated and analyzed by RNA sequencing. Read counts for individual RNAs were normalized to the total read counts of each RNA class. LPS-or VitD3-induced fold changes and corresponding p values were calculated relative to the control condition for cellular and EV-associated RNA. a Fold changes in EV-RNAs were plotted against the fold changes in cellular RNA in response to LPS (top graphs) and VitD3 (lower graphs) for miRNA (left), miscRNA (middle), and snoRNA (right). Dashed lines indicate log2FC larger or smaller than 1, so all data points beyond these lines are differentially expressed with log2FC > 1. Colored dots indicate transcripts that changed with non-adjusted p values < 0.05 in cells (blue), in EV (red), or in both cells and EV (purple), grey dots indicate unchanged transcripts. Dot size represents the normalized abundance (logCPM) of individual RNAs, mean of n = 3 experiments. b-d Transcripts showing significant changes in stimulated versus control EV (p value < 0.05) were selected, and fold changes in EV versus corresponding cells were compared. Indicated are percentages of transcripts categorized as 'reflective' (p value < 0.05 in EV and cells), 'EV-responsive' (p value < 0.05 in EV but not in cells), and 'reciprocal' (p value < 0.05 in cells and EV, but with opposite fold changes). e-g Analogous to b-d, transcripts showing significant changes in cells (p value < 0.05) were selected, and fold changes in cells were compared with those in EV. Indicated are percentages of transcripts categorized as 'reflective', 'reciprocal', cell-responsive (p value < 0.05 in cells but not in EV), or 'not in EV' (transcripts only found in cells). h RT-qPCR validation of six genes that were found to be reflective (left panels) or EV-responsive (right panels) on EV-associated RNA and cellular RNA from control, LPS-or VitD3-treated DC. Fold changes in EV and in cells were calculated relative to the control condition in cells or EV. N = 4, one-way ANOVA, *p < 0.05 ◂ 1 3 data implicate that, although the presence of specific EV-RNAs may correlate with disease, the EV transcriptome is not necessarily predictive for RNA levels in the EV-producing cell.
We observed two classes of RNA, i.e., snRNA and tRNA, of which the EV-associated levels remained stable after differential stimulation of DC. These RNA classes have been commonly detected in small RNAseq studies of EV from multiple cell types and body fluids, but factors impacting their incorporation into EV had until now not been investigated. Although these data need confirmation in different cell types with various external stimuli, our data may be a first indication that tRNAs and snRNAs are constitutive components of EV. Speculatively, these RNAs may be required for basal scaffolding functions in EV biogenesis, or may aid the functional transfer of gene-regulatory RNAs. Based on their stable association with EV, these types of RNAs also have potential to be used as EV-reference RNAs for RT-qPCR normalization, since comparing the EV transcriptome of differentially stimulated cells is complicated by potential differences in the number of EV, RNA yield, and RNA composition. Importantly, the use of well-known stably expressed cellular RNAs is not necessarily valid for normalization of EV qPCR data due to potential differences in RNA shuttling between tested conditions. Therefore, normalization of RT-qPCR data to reference genes which are stable in EV is a reliable way to quantitate differences in individual RNAs. Here, we show how RNA sequencing and RT-qPCR validation can be used to identify candidate EV-reference genes which are constitutively expressed between differently stimulated cells.
In conclusion, this study indicates that external stimuli imposed on cells can lead to changes in the levels of some EV-associated RNA classes, while leaving others unchanged. In our DC-model, two different immunerelated cell stimuli led to disparate changes in the type and quantity of miRNAs, snoRNAs, and Y-RNAs in EV. Only part of the changes in RNA contents of EV reflected changes in the cellular RNA, which urges caution in interpreting EV as snapshots of cells. These findings show that differences in multiple structurally distinct RNA classes shape the EV transcriptome. These RNA types may contribute to the functional properties of EV in cellular communication, and may be further explored as candidate EV-RNA-based indicators for the immune status of the EV-producing cell.