Characterization and Transcriptome Analysis of Exosomal and Nonexosomal RNAs in Bovine Adipocytes

Exosomes are endosome-derived extracellular vesicles that allow intercellular communication. However, the biological significance of adipocyte exosomal RNAs remains unclear. To determine the role of RNAs from bovine adipocytes and exosomes in bovine adipogenesis, exosomal and nonexosomal RNAs were extracted from three bovine primary white adipocyte samples and then profiles were generated using DNBSEQ/BGISEQ-500 technology. The RNAome of adipocytes consisted of 12,082 mRNAs, 8589 lncRNAs, and 378 miRNAs for a higher complexity that that detected in exosomes, with 1083 mRNAs, 105 lncRNAs, and 48 miRNAs. Exosomal miRNA-mRNA and lncRNA–miRNA–mRNA networks were constructed and enrichment analysis was performed to predict functional roles and regulatory mechanisms. Our study provides the first characterization of RNAs from bovine adipocyte and exosomes. The findings reveal that some RNAs are specifically packaged in adipocyte-derived exosomes, potentially enabling crosstalk between adipocytes and/or other cells that is mediated by exosomes. Our results greatly expand our understanding of exosomal RNAs from bovine adipocytes, and provide a reference for future functional investigations of adipocyte exosomal RNAs under normal physiological conditions.


Introduction
In mammals, adipose tissues serve as the largest energy storage organs, playing crucial roles in whole body energy homeostasis by secreting multiple active substances, such as adipokines [1]. Adipose tissues are composed of mostly mature adipocytes, as well as stromal cells such as preadipocytes, fibroblasts, endothelial cells, and immune cells [2]. In general, there are two main kinds of adipocytes: white adipocytes and brown adipocytes. White adipocytes can recycle and release fatty acids through lipogenesis and lipolysis, respectively, and brown adipocytes consume glucose and lipids to modulate thermal homeostasis [3]. Alterations in the mass of white adipose tissue (WAT) result from changes in white adipocyte number and/or size, in processes known as hyperplasia and hypertrophy, respectively. Adipocyte hyperplasia is related to a complex crosstalk between the proliferation and differentiation of preadipocytes, and this crosstalk is tightly regulated by both cell cycle regulators and differentiating factors. The synthesis or uptake of fatty acids leads to the accumulation of lipids in adipocytes, resulting in hypertrophy [4]. Beyond energy metabolism, the function of WAT as a dynamic secretory organ with pleiotropic functions was recently demonstrated, as diverse factors released by adipocytes can directly act on cells (autocrine) or interact with neighboring (paracrine) or distant (endocrine) cells to affect a variety of biological and pathological processes [5].
Accumulating evidence suggests that cells can secret exosomes, a kind of nanoscale membrane vesicles (<100 nm) that protect contents such as nucleic acids, proteins, and lipids from degradation and facilitate intercellular transmission, altering behaviors of recipient cells [6,7]. Exosomes are formed by inward budding of the plasma membrane, and recipient cells acquire exosomes through a range of endocytic pathways, including phagocytosis, clathrin-dependent endocytosis, and caveolin-mediated uptake [8,9]. Exosome content is highly heterogenous, depending on the origin and physiopathologic conditions of donor cells [10]. There is growing interest in the discovery and characterization of exosomal RNAs, including coding RNAs and noncoding RNA species such as miRNAs and lncRNAs, as some of these RNAs have specific functions and are selectively enclosed by exosomes in donor cells [11,12]. Adipocytes can also secrete exosomes, and there is growing interest in studying the function of adipocyte exosomal RNAs, however, the distribution characteristics and comprehensive expression profile of RNAs from adipocyte-derived exosomes remain poorly understood [13].
In this study, we systematically characterized exosomal and nonexosomal RNA profiles during the proliferation of bovine primary white adipocytes using deep sequencing technology. Based on the obtained RNA profiles, the potential sources of exosomal RNAs were determined and the potential biological functions were predicted using bioinformatics analysis. These findings provide insight into the physiological functions of bovine adipocyte-derived exosomal RNAs, and these results deepen our understanding of the biological significance of these RNAs in adipogenesis in mammals.

Identification of Adipocyte Exosomes in Cell Culture Medium
To profile exosomal RNAs, we extracted exosomes from cultured medium of bovine adipocytes and prepared RNA-seq libraries for high-throughput sequencing (Figure 1a). TEM examination of the samples revealed exosome pellets that were spherical in structures, with diameters of nearly 100 nm ( Figure 1b). NTA measurement indicated that most particles were between 60 and 80 nm in size (Figure 1c), and the particle concentration at a 5-fold dilution was approximately 8 × 10 9 particles/mL (Figure 1d), consistent with typical exosome morphology. As shown in Supplementary Materials Figure S1, Western blot analysis determined the expression of exosome surface biomarkers (CD63 and Alix).

Overview of RNA Sequencing Data
DNBSEQ/BGISEQ-500 technology was applied to systematically identify exosomal and nonexosomal RNAs from bovine adipocyte-derived exosomes and bovine adipocytes, respectively. A total of 12 sequencing libraries were obtained (including six small RNA libraries). High-throughput sequencing produced a mean of 110.8 million reads per sample, and small RNA sequencing showed an average of about 21.6 million tags per library. The q20 values of clean data for these two libraries were >97% and >99%, respectively, indicating the high quality of sequencing results (Tables S1-S3). Except for an outlier (adipo-exo3), the fractions of adipocyte samples were more similar to each other than to the samples of the adipo-exosome fractions and vice-versa (Figure 2a,b). Principle component analysis (PCA) was used to analyze the distribution and uniformity of the independent samples. Similar RNA profiles were obtained for adipocytes and adipo-exosomes (except for adipo-exo3) (Figure 2c). A boxplot was used to evaluate the overall quality of the data, including the distribution of RNA expression levels for each sample and the degree of dispersion of the data distribution ( Figure 2d). From these data, we see that there is a relatively low abundance of small RNAs in the exosomes. Additionally, our results showed 20.3% annotated miRNAs in adipocytes and 0.01% annotated miRNAs in adipo-exosomes (Figure 2e, Table S4). As shown in Figure 2e, the percentages of snRNAs and snoRNAs were low in both adipocytes and adipo-exosomes, with abundant rRNAs and tRNAs in adipo-exosomes. We also analyzed the length and frequency distribution of total small RNAs and found that most RNAs in adipocytes were 20-25 nt in length, while longer RNAs were found in adipo-exosomes, with a peak in length ≥ 30 ( Figure 2f, Table S5). These results suggest there is a mechanism to control the sort of RNA put into exosomes.

RNA Profiling in Adipocytes and Released Exosomes
In our analysis, adipo-exo3 was removed as an outlier (as described above) and we focused on RNAs such as mRNAs, lncRNAs, and miRNAs. We compared the RNA sets obtained from adipocytes and from released exosomes, and observed 1236 RNAs (Figure 3a) that were commonly expressed in adipocytes and released exosomes, including 1083 mRNAs, 105 lncRNAs, and 48 miRNAs (Figure 3b). A total of 21,049 RNAs (12,082 mRNAs, 8589 lncRNAs, and 378 miRNAs) were only detected in adipocytes and were not be expressed in released exosomes, indicating adipocytes exported only a subpopulation of RNAs into exosomes and that the process of RNA packing into exosomes is selective (Figure 3a,b). In addition, five mRNAs (EFL1, CPS1, PCM1, SH3D19, and XM_015468604.1) and three lncRNAs (BGIR9913_49345, BGIR9913_54344, and URS0000B2F7C9) were detected in exosomes irrespective of cellular origin, intimating potential different mechanisms in sorting of RNA cargoes. As shown in Figure 3c, the 12 most abundant mRNAs, lncRNAs, and miRNAs were separately filtered out from adipocytes and released exosomes. Comparison of these filtered RNAs showed adipocytes and adipo-exosomes exhibited similar profiles of abundant RNAs, with some of the same mRNAs (TMSB4X, RPS27, RPS12, TPT1, and S100A4), lncRNAs (URS0000B2207B, URS0000B2814A, XR_001499553.1, BGIR9913_54371, and XR_001499552.1), and miRNAs (bta-miR-127, bta-miR-24-3p, bta-let-7i, and bta-let-7a-5p) (Figure 3c). Applying the criteria of fold change ≥2, differential expressed RNAs between parental adipocytes and exosomes were identified, including a total of 498 mRNAs (386 upregulated and 112 downregulated), 68 lncRNAs (62 upregulated and six downregulated), and 46 miRNAs (all downregulated). Volcano plots and hierarchical clustering consistently showed variation of lncRNA, mRNA, and miRNA expression between parental adipocytes and exosomes (Figure 3d,e). All of these data revealed altered RNA profiling between adipocytes and released exosomes and indicated the selective assembly of RNAs in exosomes secreted from adipocytes rather than a random process.

Functional Enrichment Analysis of Exosomal mRNAs
Next, the top 500 abundant exosomal mRNAs were filtered out and functional enrichment was assessed. Gene ontologies (GO) analysis revealed enrichment of 49 Biological Process (BP), 84 Cellular Component (CC), and 57 Molecular Function (MF) terms with Q < 0.05 (Tables S6-S8). For the most significantly enriched components, BP terms included translation, cytoplasmic translation, and protein folding; CC terms of ribosome, cytosolic small ribosomal subunit, and cytosolic large ribosomal subunit; MF terms included structural constituents of ribosome, nucleotide binding, and protein binding ( Figure 4a). Applying the criteria of p-value < 0.05, analysis revealed 34 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched (Table S9), with many pathways related to the ribosome, biosynthesis of amino acids, and protein processing in endoplasmic reticulum ( Figure 4b). Further, reactome pathway enrichment analysis showed significantly affected pathways included formation of a pool of free 40S subunits, eukaryotic translation termination, and peptide chain elongation ( Figure 4c, Table S10). Overall, these results suggested that the selected exosomal mRNAs function in the process of protein synthesis.

Construction and Analysis of the miRNA-mRNA Network
After processing, 11 miRNAs and 328 mRNAs were identified and used to assemble the miRNA-mRNA network ( Figure 5a). In this analysis, each node represents an RNA, and each edge represents an interaction between two RNAs. The node degree typifies the number of edges linked to a node, where the higher the degree, the bigger the size, indicating biological functions. GO analysis revealed 633 BP, 45 CC, and 115 MF terms with Q < 0.05 (Tables S11-S13). By this analysis the most significant terms were BP terms including negative regulation of apoptotic process, positive regulation of gene expression, and protein phosphorylation; CC terms of cytoplasm, nuclear chromatin, and nucleus; and MF terms of protein kinase activity, protein kinase binding, and RNA polymerase II proximal promoter sequence-specific DNA binding (Figure 5b). Applying the criteria of Q-value < 0.05, KEGG analysis revealed enrichment of 88 pathways (Table S14), with many of these pathways related to the FoxO signaling pathway, the PI3K-Akt signaling pathway, and the cell cycle (Figure 5c). These results suggested that selective exosomal miRNAs act in processes of cell proliferation and apoptosis.

Coexpression of the lncRNA-miRNA-mRNA Network and Functional Prediction
As shown in Figure 6a, 24 lncRNA, 15 miRNA, and 98 mRNA composed the lncRNA-miRNA-mRNA network (Figure 6a). The results of the GO analysis revealed 493 BP, 26 CC, and 89 MF terms with Q < 0.05 (Tables S15-S17). The most significant enrichment was for BP terms of regulation of gene expression, positive regulation of apoptotic process, and regulation of cell proliferation; CC terms including cytosol, nucleus, and nuclear chromatin; MF terms of transcription factor binding, nucleotide binding, and protein kinase activity (Figure 6b). Applying the criteria of Q-value < 0.05, 50 KEGG pathways were enriched (Table S18), with many of these pathways related to the FoxO signaling pathway, cellular senescence, and the p53 signaling pathway (Figure 5c). These data show that bovine adipocytes released exosomal lncRNAs that are predicted to regulate cell-based proliferation and apoptosis.

Discussion
Since the initial identification of leptin and adiponectin, WAT has been recognized as a dynamic secretory organ with pleiotropic functions beyond energy metabolism [14]. WAT includes a considerable variety of cell types, including preadipocytes, fibroblasts, immune cells, and endothelial cells [15][16][17][18]. Cells secrete diverse factors and adjacent cellular communication is mediated by gap junctions or tunneling nanotubes [19]. The recent characterization of exosomes and their carried contents has led to an extension of our understanding of intercellular communication. Exosomes provide a natural carrier system that can transfer carried contents such as proteins, nucleic acids, and lipids between donor and receiver cells, acting to protect loaded cargoes from being eliminated by the mononuclear phagocyte system and providing a flexible mechanism of intercellular transfer [20]. Several reports have characterized the potential RNA contents of exosomes, and exosomal lncRNA and miRNA have demonstrated functions in recipient cells. However, only preliminary characterization and analysis of exosomes and carried RNAs in adipocytes have been performed. In this work, we isolated and verified bovine white adipocyte-derived exosomes; comprehensively profiled the mRNA, lncRNA, and miRNA transcriptomes of exosomes using RNA-seq; ultimately evaluated the molecular relatedness of these RNAs and analyzed their putative functions.
Our data provided the first insight into differences between bovine adipocytes and released exosomes, and results supported the hypothesis of directed RNA export in exosomes. Both Pearson correlation analysis and PCA showed highly correlated RNA profiles between adipocytes and adipo-exosomes, except for one adipo-exo3, so this one was removed as an outlier in our subsequent analysis. The RNA sequencing revealed generally higher expression of exosomal mRNAs and lncRNAs than that of cellular mRNAs and lncRNAs, and lower miRNA expression levels in exosomes compared to that in adipocytes, indicating differential sorting of various RNAs into exosomes. Several recent studies have reported mechanisms involved in the selective targeting of RNAs into exosomes. For example, some motifs enriched in mRNA 3 -UTR may potentially act as cis-acting elements to promote targeting of mRNAs into exosomes [21], and targeting lncRNA to exosomes occurs via a specific protein hnRNPA2B1, which characteristically binds to the 5 sequence of lncARSR to guide lncRNA sorting into exosomes [22]. The packaging of miRNAs into exosomes can be mediated by the neural sphingomyelinase 2 (nSMase2)-dependent pathway [23] and 3 end uridylation [24], as well as by RNA-binding proteins such as Ago2 [25], heterogeneous nuclear ribonucleoproteins (hnRNP) family protein [26], and Y-box protein 1 (YBX1) [27]. Specific lncRNA-RBP complexes can capture special miRNAs and sort them into exosomes [28]. As the most common fragments ranging in size from 19 to 25 nt, miRNAs make up the largest fraction of RNA in cells, however, only 0.01% miRNAs were detected in exosomes. Interestingly, abundant rRNAs and tRNAs were found in adipocytes and released exosomes. Similar to previous reports of exosomal RNA enrichment [29,30], we found that exosomes were enriched with 5S ribosomal RNA in the absence of visible 18S and 28S rRNA peaks on exosomal RNA profiles ( Figure S2). Another abundant RNA species, tRNAs represented over 50% of the total small RNAs in adipose-derived exosomes. Although it is unknown whether these exosomal rRNAs and tRNAs have a specific function, some limited evidence has suggested the presence of fragmented rRNA and tRNA in exosome, which may represent novel RNA types with potential regulatory functions [30,31]. The enrichment of small RNAs also exhibited a similar size distribution as that previously reported in bovine sera and exosomes with a ≥30 peak observed in adipo-exosomes; this peak may correspond to processed small RNAs acting as a novel form of signaling molecules [32]. However, further investigation is required to determine whether and how these different classes of exosomal RNA influence recipient cells.
In this study, a total of 1236 RNAs were normally expressed in adipocytes and released exosomes, including 1083 mRNAs, 105 lncRNAs, and 48 miRNAs. Although fewer in number, these exosomal RNA profiles were closely correlated with those obtained previously for native cells. This is consistent with a wide array of studies suggesting that cell-derived exosomal RNA profiles closely reflect RNA profiles in parental cells, and these exosomal RNAs may transfer intercellular signals to affect critical biological functions [33,34]. First detected in murine and human mast cell lines (MC/9 and HMC-1), exosomal mRNAs have been described as important regulators of cell behavior, with approximately 1300 mRNAs detected from MC/9 exosomes by microarray analysis and associated with ontologies such as RNA post-transcriptional modification and protein synthesis. Interestingly, the transfer of murine exosomal mRNA to the HMC-1 generates the output of new murine proteins in the HMC-1, indicating that exosomal mRNAs can be translated into proteins in recipient cells [35]. We observed enrichment in exosomes for mRNAs that encode ribosomal proteins in the small (RPS) or large (RPL) subunits, such as RPS2, RPS12, RPS27, RPL18, RPL23, and RPL32, and functional annotations of exosomal mRNAs have been consistently linked to the biosynthesis of amino acids and protein processing and modification [35,36]. Exosomal mRNAs released from adipocytes can be transmitted to recipient cells, where they can facilitate the protein expression programs of recipient cells. Jenjaroenpun et al. also hypothesized that exosomal mRNAs can be translated into proteins that support ribosomal function after being delivered to recipient cells, facilitating the efficacious translation of other exosomal mRNAs in a cellular domain where exosomal content is released [30].
Our results showed only a small number of miRNAs in adipocyte exosomes, but the few miRNAs detected may be important in adipocyte biology. Specifically, expression of miR-143, miR-15b, and let-7 subtypes was expressed in both adipocytes and released exosomes. Ectopic expression of miR-143 in preadipocytes accelerates the rate of adipocyte formation and acts in adipocyte differentiation by targeting Erk5 [37]. In porcine preadipocytes, the inhibition of FOXO1 expression by miR-15b indicates a positive effect on adipogenesis [38]. The let-7 miRNA modulates the larva to adult transition in caenorhabditis elegans, and a study on 3T3-L1 cells supports a role for let-7 in the transition to terminal differentiation during adipogenesis [39,40]. With high conservation across species and wide regulatory roles in target gene expression, these miRNAs may directly influence cell functions once released into other adipocytes. Interestingly, there is no report on the effects on adipogenesis by these abundant exosomal miRNAs, but there is potential for significant regulation of surrounding cells. As major exosomal signals and effectors, a broad range of adipocyte-derived miRNAs can transfer into recipient cells via exosomes, resulting in the suppression of mRNA targets. For example, miR-27a can repress adipogenesis by targeting PPARγ, and recent studies demonstrated that adipocyte-derived exosomal miR-27a could be assimilated by myocytes triggering insulin resistance in skeletal muscle through repression of PPARγ [41]. Pan et al. identified adipocyte-secreted miR-34a as a key paracrine mediator that suppresses M2 polarization by repressing KLF4 in adipose-resident macrophages [13]. Another miRNA, miR-23a/b is delivered from adipocytes to cancer cells, and miR-23a/b was shown to act to modulate HCC cells by repressing the expression of VHL/HIF-1α [42]. The functional exosomal miRNAs described above were all found in bovine adipocyte-released exosomes, which increase the credibility of our sequencing results to a certain degree. Similarly, the abundant exosomal miR-451 detected in this study has also been reported as a highly exported exosomal miRNA in many different cell types, where it has been shown to be processed by the Ago2 protein rather than RNase III Dicer into mature miRNA duplexes [43,44]. To identify functions of abundant exosomal miRNAs, GO-based target prediction was performed and revealed potential roles for these miRNAs in apoptotic process, regulation of cell proliferation, and protein phosphorylation, with significant enrichment of pathways related to FOX and PI3K-AKT3 signaling. These predictions suggest possible activities of these exosomal miRNAs in cell proliferation and apoptosis, however detection in exosomes does not necessarily mean that all these miRNAs are transferred into other cells, so this requires further investigation.
As mRNA-like transcripts over 200 nucleotides length with abundant binding sites for miRNA, lncRNAs competitively interact with miRNA to repress mRNA targets expression, in a mechanism known as competing endogenous RNA (ceRNA) [45]. With low conservation across species, almost all of our identified lncRNAs have not been reported in previous studies. Accordingly, we constructed an exosomal lncRNA-miRNA-mRNA network for exosomal lncRNAs annotation. The ceRNA network has been shown to regulate a variety of cellular processes in adipocytes. For example, lncRNA ADNCR can inhibit adipocyte differentiation via sponging miR-204, indirectly augmenting gene expression of SIRT1, which can inhibit adipocyte formation [46]; GAS5 can compete for miR-18a/CTGF to suppress adipogenic differentiation of mesenchymal stromal cells (MSCs) [47]; RP11-142A22.4/miR-587/Wnt5β axis is a potential pathway to promote adipogenesis through the ceRNA mechanism [48]. However, there are few ceRNA studies about exosomal lncRNAs. The network analysis results indicated that BGIR9913_50763, BGIR9913_55051, BGIR9913_51021, and BGIR9913_52429 are crucial components of this network, with links to the greatest number of miRNAs. Furthermore, functional enrichment analysis revealed significant enrichment of these exosomal lncRNAs in cell proliferation and apoptosis-related biological processes, including the FOX and p53 signaling pathways, suggesting their importance in cell growth and development.
As a natural cell product, exosomes transport functional RNAs between cells, avoiding their degradation by phagocytes or lysosomes [49]. Hence, exosomes are logically considered as key contributors to mechanisms that adipocytes can employ to influence other cell systems. In this study, we investigated the RNA transcriptome of bovine adipocytes and their released exosomes, using deep sequencing technology, expanding our understanding of the physiological functions and mechanisms of adipocyte exosomal RNAs. To our knowledge, this is the first systemic study of adipocyte exosomal RNAs in bovine, and the results of this work should provide a reference for future functional studies on adipocyte exosomes under normal physiological conditions. However, there are still some limitations of this work that should be acknowledged. Given the current limitations of exosome separation techniques, it is difficult to completely discriminate exosomes and specific subtypes from other vesicle types, resulting in the heterogeneity of RNAs in exosomes [9]. Therefore, our exosomal RNAs in the libraries may be merely subsets of all RNAs isolated in this experiment. Overall, improvements in the stringent purification of exosomes are required. Similarly, it is difficult to separate high-purity adipocytes from WAT, raising the possibility that exosomes from fibroblasts, endothelial cells, or immune cells may contaminate our exosome isolation. In addition, reads from RNA-seq are calculated based on signal intensities, which are not always accurate due to technical limitations, so biological verification should be performed to validate the results of sequencing. Recently, coculture systems were used to reproduce cellular interactions, with exosomal RNAs produced by overexpression of exosomes [50,51]. However, it remains unclear whether these overexpressed exosomes can produce functional RNAs that can be transferred in native biological settings, and the slight induction of exosomal RNAs observed in target cells may not be sufficient to produce significant effects. Furthermore, with some common targets RNAs in the exosomes may collaborate to induce target cell phenotypes. These issues will need to be addressed in future work.

Ethics Statement, Sample Collection, and Cell Culture
All animal experiments were conducted according to the Regulation for the Administration of Affairs Concerning Experimental Animals (Shaanxi, China, 2017), and our study was approved by the Northwest A&F University Institutional Animal Care and Use Committee (November 2017). Adipose tissues were acquired from Qinchuan cattle at the embryonic stage (90 days), and these samples were provided by Shannxi Kingbull Animal Husbandry (Baoji, China). Adipocytes were extracted from bovine adipose tissues as previously detailed [52,53], and were then seeded into 100-mm plates and cultured at 37 • C with 5% CO 2 until they reached 90% confluence using DMEM medium (20% exosomes-depleted fetal bovine serum with 1% penicillin-streptomycin). Three biological replicates were performed for adipocyte culture.

Exosome Isolation
The conditioned medium was collected in a fresh 50 mL tube and centrifuged at 800 rpm for 10 min followed by 12,000 rpm for 30 min at 4 • C to remove debris. After centrifugation, the supernatant was filtered using a 0.22-µm filter (Millipore, Burlington, MA, USA) to eliminate large membrane vesicles, and then was further centrifuged at 100,000 rpm for 70 min to deposit exosomes. For subsequent experimental verification, exosomes were obtained using the ExoQuick solution (SBI, Fremont, CA, USA) performed according to the manufacturer's instructions, with exosomes resuspended in elution buffer and stored at 80 • C before further analysis.

Identification of Bovine Adipocyte Exosomes
To evaluate the exosome characteristics, nanoparticle tracking analysis (NTA) was performed using a Flow NanoAnalyzer (NanoFCM, Fujian, China) according to the manufacturer's instructions. In addition, the morphology of exosome pellets was observed by transmission electron microscopy (TEM) HT-7800 (HITACHI, Tokyo, Japan), and exosomal surface markers (CD63 and Alix) were detected via Western blotting analysis. All samples were evaluated in triplicate.

RNA Extractions and Quantitation
Total RNA was extracted from bovine adipocytes and exosomes using the Trizol reagent (TAKARA, Shiga, Japan) and RA806TC-1 (System Biosciences, Palo Alto, CA, USA) following manufacturer's protocols. The RNA quantity and concentration were assessed using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA), and the integrity was determined by agarose gel electrophoresis. Total RNA samples were stored at −80 • C before further analysis.

RNA Library Construction and Deep Sequencing
The construction and deep sequencing of RNA libraries were accomplished with the assistance of Shenzhen BGI Co., Ltd. Briefly, RNA samples from bovine adipocytes (n = 3) and exosomes (n = 3) were ligated with 5 and 3 adapters, followed by cDNA synthesis using adaptor-specific primers. After PCR amplification, the RNA libraries were gel purified and the library quality was assessed. Then, sequencing was performed on the DNBSEQ/BGISEQ-500 platform (BGI, Shenzhen, China).

RNA Sequence Analysis
Clean reads were saved in FASTQ format, and low quality reads were removed from raw data using SOAPnuke (https://github.com/BGI-flexlab/SOAPnuke). Valid reads were aligned against both the bovine miRNA sequences (Release 21) and the reference bovine genome (GCF_000003205.7_Btau_5.0.1) to annotate the known RNAs using HISAT2 (http://www.ccb.jhu.edu/software/hisat) and Bowtie2 was used for miRNA. RSEM (http://deweylab.biostat.wisc.edu/rsem) was used to estimate RNA abundances, and fragments per kilobase of transcript per million mapped reads (FPKM) were determined for RNA transcript expression levels in this study. Here, FPKM > 0 was used as a value of cut-off limit for RNA detection across all samples.

Functional Enrichment Analysis
For the functional exploration of exosomal mRNAs, the top 500 abundant exosomal mRNAs were screened, and functional annotations including GO analysis, KEGG pathway analysis, and reactome enrichment analysis were performed using Phyper function in R software. Q-value < 0.05 (p-value < 0.05 for KEGG pathway analysis) was used as a threshold of significance.

Construction and Analysis of the Network of Exosomal RNAs
In our study, the 12 most abundant miRNA were chosen from 48 exosomal miRNA, and their target genes were further predicted using miRTarBase, an experimentally confirmed miRNA-mRNA interactions database (Chou et al., 2018). These acquired miRNA-mRNA pairs were visualized using Cytoscape software (http://www.cytoscape.org/) and GO and KEGG pathway analysis were performed for exosomal miRNA function annotation. Previous work has shown that lncRNAs can regulate gene expression by acting as miRNA sponges as part of the competing endogenous RNAs network. TargetScan (http://www.targetscan.org/), miRanda (http://www.microrna.org/microrna/home.do), and RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid) were used to predict lncRNA-miRNA target relationships. The target miRNAs of exosomal lncRNAs were further screened as MRE < −36 standard, and miRNA-mRNA pairs were assessed by bioinformatics analysis using miRTarBase program [54]. Ultimately, the lncRNA-miRNA-mRNA network was reconstructed and visualized by Cytoscape software (Cytoscape 3), and GO and KEGG pathway analyses were performed again for exosomal lncRNA function annotation.

Conclusions
Here, for the first time, we isolated and characterized exosomes derived from bovine adipocytes, and identified mRNAs, lncRNAs, and miRNAs with the potential to regulate recipient cell phenotype and modulate multiple cellular pathways. We detected a notable enrichment of RNAs in exosomes compared to parental cells, and the differences in RNAs suggest a complex sorting mechanism of RNAs into exosomes. The results of this work provide a better understanding of exosomal RNAs in bovine adipocytes and should facilitate further study of bovine adipogenesis and exosomal RNAs.
Supplementary Materials: Supplementary Materials can be found at http://www.mdpi.com/1422-0067/21/23/ 9313/s1. Figure S1: Western blot analysis of exosomes derived from bovine adipocytes. Adipocytes lysates and exosomes were blotted for CD63 and ALIX. Figure S2: Bioanalyzer profiles of the total RNAs of bovine adipocytes and released exosomes. Table S1: Quality statistics of filtered reads. Table S2: Quality statistics of filtered reads (small RNA). Table S3: Alignment clean reads to the reference genes. Table S4: smallRNA Classification.  Table S5: smallRNA length distribution. Table S6: BP result of the exosomal mRNAs. Table S7: CC result of the exosomal mRNAs. Table S8: MF result of the exosomal mRNAs. Table S9: KEGG pathway result of the exosomal mRNAs. Table S10: Reactome Enrichment result of the exosomal mRNAs. Table S11: BP result of the miRNA-mRNA network. Table S12: CC result of the miRNA-mRNA network. Table S13: MF result of the miRNA-mRNA network. Table S14: KEGG pathway result of the miRNA-mRNA network. Table S15: BP result of the lncRNA-miRNA-mRNA network. Table S16: CC result of the the lncRNA-miRNA-mRNA network. Table S17: MF result of the lncRNA-miRNA-mRNA network. Table S18: