Isolation and identification of exosomes from anterior pituitary of Duroc swine
Exosomes were isolated from Duroc swine anterior pituitary (Additional Fig.1). We detected the purified vesicles using transmission electron microscopy which showed that their size and morphology in cup-shaped (Fig.1a), typical of that of exosomes. Then we used the Zetasizer to analyze their size distribution and found that the vesicles average size is about 92nm (Fig.1b). Moreover, exosomes were confirmed by Western Blot with positive common surface markers CD9 and CD63 (Fig.1c).
Overview of small RNA deep sequencing data in exosomes and analysis
In order to explore the non-coding RNA expression profiles of the exosomes from anterior pituitary of Duroc swine breed, we used RNA-seq analyses to characterize the non-coding RNA from three normal anteriors of 60-day-old Duroc swine. We obtained 12778982 (Exo_1), 15668033 (Exo_2) and 1535301139 (Exo_3) clean reads that were screened from sRNA for subsequent analysis after quality evaluation (Additional file 1: Table S1). Meanwhile, the length distribution of the obtained total sRNA fragment was statistically analyzed (Additional Fig.2). In general, sRNAs range from 18 to 35nt in length and the majority of the miRNA reads were about 22 nt. A total of 416 miRNAs were obtained from samples, 343 of which are known miRNAs and 73 are newly predicted miRNAs (Additional file 2: Table S2). Of these known miRNAs, 61 miRNAs were highly expressed (1,000 < average signals ≤ 10,000), and, in particular, 46 miRNAs were extremely highly expressed in exosomes (average signals ≥ 10,000). To further characterize the regulatory roles of miRNAs in exosomes from anterior pituitary, miRNA target prediction, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation analyses were performed in our study. A total of 185183 target genes for the 416 miRNAs were predicted. Our GO annotation indicated that the predicted target genes were significantly enriched in intracellular signal transduction, phosphorylation, catabolic process, developmental process, the component of cytoskeletal part, binding, protein binding and nucleotide binding (Fig.2a and Additional file 3:Table S3). The significantly enriched KEGG pathways mainly included PI3K−Akt signaling pathway and Calcium signaling pathway (Fig.2b and Additional file 4: Table S4). Previous studies documented that the PI3K/Akt singling pathway participate in regulating growth hormone secretion[84] and Calcium signaling pathway involved in controlling excitability of anterior pituitary cell[85]. These findings suggest that miRNA in exosomes of anterior pituitary could regulate intracellular signal transduction, catabolism and development.
Overview of lncRNA deep sequencing data in exosomes and bioinformatics analysis
LncRNA is a class of RNA molecules with transcript length over 200nt and do not encode proteins. According to its characteristics, we set filter criteria and counted the number of transcripts screened per step (Additional Fig.3a). For lncRNAs prediction, CPC and CNCI were used for potential coding ability detection, and PFAM, a protein database, was used for protein annotation information analysis were used for potential coding ability detection (Additional Fig.3b). Resultantly, 15545 novel lncRNAs and 687 annotated lncRNAs (Additional file 5:Table S5) were identified respectively. We performed statistics on different types of lncRNA mainly for lincRNA, anti-sense_lncRNA and intronic_lncRNA. The results showed that the percentage of intronic_lncRNA was the highest (Fig.3a). The structure and sequence conservation of lncRNA and mRNA were also compared and analyzed. We found that lncRNAs were with shorter length in transcript (Additional Fig.3c) and their genes tend to contain fewer exons (Fig.3b). Most of the mRNA had longer open reading frames than lncRNA (Additional Fig.3d). The transcript expression level of lncRNA was higher than miRNA (Additional Fig. 3e) and we got same perception by comparing the FPKM of exosomes from the different samples (Additional Fig. 3f). We investigated the possible functions of the lncRNAs by searching for protein-coding genes 100 kb upstream and downstream of all identified lncRNAs to predict the potential cis-regulatory targets of lncRNAs. A total of 57087 protein-coding genes were closest to the 9524 lncRNAs. A number of lncRNAs were detected co-expression with pituitary-specific genes including growth hormone 1 (GH1), growth hormone releasing hormone receptor (GHRHR), prolactin releasing hormone receptor (PRLHR), follicle stimulating hormone subunit beta (FSHB), luteinizing hormone subunit beta (LHB) (Fig.3c). Otherwise, some lncRNAs could co-expressed with genes involved in exosome marker protein, protein transport and exosome docking such as CD63, CD81, TSG101, Rab27A, Rab27B, UBL3. Our GO annotation indicated that the predicted target genes of lncRNA were significantly enriched in cellular component biogenesis, organelle organization, RNA biosynthetic process, the cellular component of nucleus and organelle part, organic cyclic compound binding, nucleotide binding and small molecule binding (Fig.3e and Additional file 6:Table S6). The KEGG pathways enriched in Ras signaling pathway, Hippo signaling pathway and Cell cycle (Fig.3f and Additional file 7:Table S7). Hippo signaling pathway could play a role in pituitary development and regulating pituitary stem cells[86]. Ras signaling pathway could regulate pituitary cell-specific gene expression[87]. These data indicate that lncRNA in exosomes of anterior pituitary could participate in the development of pituitary and RNA biosynthetic process.
Overview of circRNA deep sequencing data in exosomes and bioinformatics analysis
After evaluated the data output quality, we obtained 494 novel circRNAs (Additional file 8:Table S8) and then counted the length distribution and the source of circRNA for all samples (Fig.4a). It showed that the length of circRNAs are mostly scattered in the area less than 10000nt and the source of circRNA mostly from intergenic area compared with exon and intron area (Fig.4b). The expression levels of all circRNAs were statistically analyzed and normalized by TPM (Fig.4c). TPM density distribution allows overall inspection of gene expression patterns in samples and the results showed there are lots of overlap which means consistency between samples[88]. We then constructed a circRNA-miRNA co-expression network based on the RNA-seq results. CircRNA could inhibit the function of miRNA by combining with miRNA[89]. So the analysis of miRNA binding sites on the identified circRNA helps further study for the function of circRNA. Then we used miRanda software to predict the miRNA binding site of the cleaved circRNA and focused on the circRNAs which were combined with highly expressed miRNAs in pituitary and exosomes from anterior pituitary. A network map was constructed containing 39 circRNAs, 8 miRNAs and 49 relationships (Fig.4d). In order to explore the potential functions of the circRNAs in exosomes from anterior pituitary, we performed GO and KEGG pathway enrichment analysis. The results showed that the enriched GO terms were mainly associated with metabolic process, cellular biosynthetic process, binding and transferase activity (Fig.4e and Additional file 9:Table S9) and the KEGG pathways were mainly enriched in the Wnt signaling pathway, regulation of actin cytoskeleton, protein processing in endoplasmic reticulum and phagosome (Fig.4f and Additional file 10:Table S10). These findings suggest that circRNA in exosomes of anterior pituitary could regulate the cellular metabolic and biosynthetic process.
Analysis of crosstalk among lncRNA-miRNA-mRNA in exosomes
Recent studies suggested that lncRNAs could function as endogenous miRNA sponges to prevent miRNA from binding to reduce the regulatory effect of miRNAs on their target mRNA[90-92]. To further analyze the crosstalk between lncRNA, miRNA and mRNA, we predicted the interaction of them and further focused on the competitive endogenous RNAs (ceRNAs) relative with pituitary function. Then make the network that 97 lncRNAs could sponge 11 miRNAs to regulate 10 pituitary-specific genes including GH1, GHRHR, PRLHR, FSHB, LHB, proopiomelanocortin (POMC), growth hormone receptor (GHR), prolactin receptor (PRLR), gonadotropin releasing hormone receptor (GNRHR), POU class 1 homeobox 1 (POU1F1) (Fig.5a). Moreover, we performed enrichment analysis of all of dates using the GO and KEGG analysis. GO analysis revealed 273 significantly enriched terms (P < 0.05, Additional file 11:Table S11) in the categories of biological process, molecular function, and cellular components and we showed a part of terms with lots of gene numbers (Fig.5b). Its annotation indicated that they participate in intracellular signal transduction, cellular component organization or biogenesis, RNA metabolic process, localization, regulation of metabolic process, binding and regulate catalytic activity which suggest that they are involved in the body’s basic biological regulation. KEGG pathway analysis demonstrated that 169 terms (Additional file 12:Table S12) were enriched and we selected the top 20 (Fig.5c) in which the MAPK signaling pathway, GnRH signaling pathway and Dopaminergic synapse were involved in the regulatory function about information transfer about pituitary. These results suggest that the crosstalk among lncRNA-miRNA-mRNA could participate in pituitary signaling process.
Analysis of crosstalk among circRNA-miRNA-mRNA in exosomes
The current studies have proved that circRNAs could role as ceRNAs to compete for miRNA-binding sites to affect the function of miRNA [93, 94]. Therefore, analysis of interactions between miRNAs and circRNAs is helpful for further study. In the constructed potential circRNA–miRNA–mRNA associations, similarly we mainly concerned the ceRNAs relative to pituitary function. The resultant network was comprised of 188 edges among 11 miRNAs, 58 circRNAs and 10 pituitary-specific genes including GH1, POMC, GHR, GHRHR, PRLR, LHB, PRLHR, FSHB, GNRHR, POU1F1 (Fig.6a). In order to learn the potential functions of the associated non-coding RNA in exosomes from anterior pituitary, we conducted GO and KEGG enrichment analysis. GO analysis revealed 265 significantly enriched terms (Additional file 13: Table S13) and we showed some terms enriched a lot of gene numbers (Fig.6b). Our GO annotation indicated that they are involved in intracellular signal transduction, cellular component organization or biogenesis, transport, protein binding, hydrolase activity and phosphotransferase activity. Our KEGG pathway analysis demonstrated that 267 terms (Additional file 14: Table S14) were enriched and we showed the top 20 (Fig.6c) in which the MAPK signaling pathway, Prolactin signaling pathway, GnRH signaling pathway, Estrogen signaling pathway, and Dopaminergic synapse are were participated in the regulation of hormone secretion in the anterior pituitary. These findings suggested that the network among circRNA-miRNA-mRNA in exosomes play an important role in pituitary endocrine functions.