Skip to main content
  • Research article
  • Open access
  • Published:

Exploration of exosomal microRNA expression profiles in pigeon ‘Milk’ during the lactation period

Abstract

Background

Pigeon crop has the unique ability to produce a nutrient rich substance termed pigeon ‘milk’ (PM), which has functional resemblance with the mammalian milk. Previous researches have demonstrated that a large number of exosomes and exosomal miRNAs exist in mammalian milk, and many of them are associated with immunity, growth and development. However, to date, little is known about the exosomes and exosomal miRNAs in PM.

Results

In this study, we isolated the exosomes from PM and used small RNA sequencing to investigate the distribution and expression profiles of exosomal miRNAs. A total of 301 mature miRNAs including 248 conserved and 53 novel miRNAs were identified in five lactation stages i.e. 1d, 5d, 10d, 15d, and 20d. From these, four top 10 conserved miRNAs (cli-miR-21-5p, cli-miR-148a-3p, cli-miR-10a-5p and cli-miR-26a-5p) were co-expressed in all five stages. We speculate that these miRNAs may have important role in the biosynthesis and metabolism of PM. Moreover, similar to the mammalian milk, a significant proportion of immune and growth-related miRNAs were also present and enriched in PM exosomes. Furthermore, we also identified 41 orthologous miRNAs group (giving rise to 81 mature miRNA) commonly shared with PM, human, bovine and porcine breast milk. Additionally, functional enrichment analysis revealed the role of exosomal miRNAs in organ development and in growth-related pathways including the MAPK, Wnt and insulin pathways.

Conclusions

To sum-up, this comprehensive analysis will contribute to a better understanding of the underlying functions and regulatory mechanisms of PM in squabs.

Background

Unlike the most avian species, the pigeons (Columba livia) are among few representative birds which have the ability to produce a unique nutritive substance called pigeon ‘milk’ (PM) for the sustenance and growth of their squabs. PM is a curd like cheesy nutritive solution and has many functional similarities with the mammalian milk in terms of delivery of nutritional and immunogenic benefits and subsequent growth and development of young ones [1]. Despite such similarities, one contrasting fascination is that PM is produced by both female and male pigeons; however, there is dearth of scientific literature addressing the underlying mechanisms of these processes. The production of PM is highly specialized process, it is well established that, alike mammals, the prolactin (lactogenic hormone) regulates the production of PM in pigeons [2]. However, histological studies on lactating pigeon crop tissue have showed that the process is structurally orthogonal to conventional mammalian lactation since the pigeon crop is non-glandular and secretory processes seem to be absent in it [3,4,5]. Gillespie et al. in their genome-wide pigeon crop transcriptome study have investigated the underlying molecular mechanism of pigeon milk production [6]. Their findings have revealed that the differential expression of cornification related proteins and de novo lipid synthesis genes in pigeon crop during lactation contribute to the highly specialized process that leads to the production of PM. PM consists of lipid enriched keratinocytes (differentiated cornified cells) which after undergoing rapid proliferation have separated from crop sac germinal epithelium [7]. A large number of studies have indicated that the PM contains highly nutritious substances such as proteins, lipids, minerals, carbohydrates and growth factors with diverse biological activities [8,9,10]. Shetty et al. have reported a stupendously rapid growth in squabs during the lactation period, which provides manifestation for the physiological and evolutionary immensity of PM [11]. In previous functional studies, PM has reportedly stimulated the growth of ovarian cells in hamster [11], caused precocial incisor eruption and eyelid opening in newborn mice [12], and promoted the growth in chickens [13]. Gillespie and colleagues have observed the expression of immune-related gene pathways and significantly increased interferon-stimulated genes in the gut associated lymphoid tissue (GALT) of PM fed chickens [1]. However, one previously conducted study reports that the squabs failed to thrive when they were fed with the nutritional replacement of PM [14]. These results suggest that there are other essential factors other than the nutritional substances in PM that influence the development of young.

MicroRNAs (miRNAs) are evolutionary conserved small non-coding RNAs (~ 22 nt), which are widely involved in the complex cellular mechanisms such as differentiation, proliferation, metabolism and immune development through the post-transcriptional gene regulatory mechanisms [15,16,17]. In recent years, a large number of miRNAs were found to be present in various biological fluids i.e. breast milk, serum, saliva, and urine [18,19,20], and these extracellular miRNAs are principally delivered by exosomes. Exosomes are nano-sized membrane vesicles of endocytic origin and have emerged as important mediators of the intercellular communication in immune system and elsewhere by transferring miRNAs [21, 22]. In our previous studies, we have reported that the immune-related exosomal miRNAs were enriched in human [22], porcine [23] and giant panda breast milk [24]. Based on findings of these studies, we speculated and suggested that the exosomal miRNAs may play a significant role in growth and development of infants. However, to date, no study has reported the presence of miRNAs in PM, therefore, the potential physiological/functional role of exosomal miRNAs in PM remains unidentified.

Here, we conjectured that the exosomal miRNAs are enriched in PM and are involved in modulation of physiology and immunity related mechanisms in developing squabs. We isolated the exosomes from PM and examined the expression profiling of exosomal miRNAs using small RNA-seq approach. Furthermore, in order to elucidate the potential role of exosomal miRNAs in PM, we also performed a comparative transcriptomic analysis of exosomal miRNAs in PM and in breast milk of other three representative mammals i.e. human, cattle and pig.

Results

Identification of small RNAs-loaded exosomes in PM

Exosomes were isolated from PM (Fig. 1a) by ultracentrifugation and were investigated under the atomic force microscope (AFM) at the nanometer scale. The PM contained a substantial amount of exosome-like vesicles (Fig. 1b, c) with an approximate width of 150 nm and a height of 8 nm (Fig. 1d). They exhibited the similar ultrastructural morphology as that of exosomes previously identified in porcine [23] and bovine breast milk [25]. Additionally, as per recommendation of the International Society for Extracellular Vesicles (ISEV) the exosomes were further confirmed by Western blot analysis. As an exosomal marker, the tetraspanin molecule CD63 was detected, while the Tubulin (a cytoskeletal protein) was absent in ultracentrifugation pellets (Fig. 1e). We further identified that these vesicles contained a substantive quantity of small RNAs shorter than 100 nt in length, but no 18S and 28S ribosomal RNA were observed (Additional file 1: Figure S1), which confirmed the presence of small RNAs in exosomes.

Fig. 1
figure 1

Identification of exosomes in PM. a The morphology of PM. b 2D AFM image of isolated exosome vesicles in PM. c 3D AFM image of isolated exosome vesicles in PM. d The line profile of AFM image for a PM exosome vesicle. X- and Y-axes are the width and height of PM exosome vesicle, respectively. e The expression of CD63 and Tubulin in the crop tissue and exosome vesicles was determined by Western blotting. It can be seen that the isolated exosome vesicles are positive for CD63 but negative for Tubulin

Identification of exosomal miRNAs in PM

To identify the exosomal miRNAs in PM across the lactation periods, we generated a total of 95,476,045 raw reads of seven small RNA libraries (1d-1, 1d-2, 1d-3, 5d, 10d, 15d, 20d). After handling with a series of stringent filters (such as removing low-quality reads, repeated sequences, and adaptor sequences), we obtained a set of 81,219,606 (85%) high-quality reads for miRNA identification. With strict mapping criteria to the pigeon reference genome, retained reads were considered as reliable miRNA candidates (Additional file 2: Table S1). The length distribution peaked at 22~ 24 nt, which is consistent with the animal miRNAs (Additional file 3: Figure S2A).

A total of 230 pre-miRNAs giving rise to 301 mature miRNAs were identified in seven libraries after mapping to the pigeon reference genome (Columba livia, ColLiv2) (Table 1). Since, the pigeon miRNA sequences are currently unavailable in miRbase 21.0, the miRNAs were categorized into two groups as per following alignment criteria: 248 miRNAs were mappable to other two known avian species (Gallus gallus and Taeniopygia guttata) and mammalian miRNAs and are regarded as ‘conserved miRNA’ and prefixed with ‘conserve-cli’. Whereas, 53 miRNAs were unable to be annotated by known avian and mammalian mature miRNAs and satisfied the folding minimum free energy cutoff of RNA hairpin structure and therefore, are called ‘novel miRNAs’ and prefixed with ‘novel-cli’ (Additional file 4: Table S2). Not surprisingly, the read counts of conserved miRNAs were more abundant compared to those of novel miRNAs, which seems to be a common phenomenon across species [26] (Additional file 3: Figure S2B, C).

Table 1 Conserved and novel miRNAs detected in exosomes of PM

Immune and growth-related microRNAs are abundant in exosomes of PM

To better elucidate the miRNA expression in PM across the lactation periods, we conducted the hierarchical clustering analysis on the annotated conserved miRNA groups. As shown in Fig. 2a, the libraries clustered together based on the lactation time. Two major branches were defined: one representing 1, 5 and 10 days, and other representing 15 and 20 days. This different clustering pattern may correspond to different kinds of miRNAs as well as the different expression patterns of the same miRNAs. Moreover, three biological replicates in 1 day were highly clustered together, which suggested the experimental reliability and highlighted the low variation in miRNAs profiles in PM across the individuals. The two major branches had 115 miRNAs overlapped, accounting for 49.78 and 87.12% in the first three stages (1d, 5d, 10d) and the last two stages (15d, 20d), respectively. Furthermore, 121 miRNAs were shared in 1d, 5d and 10d, only a few miRNAs were uniquely expressed in 5d and 10d, which indicates the time-dependent decreasing tendency in expressed miRNAs species (Fig. 2b). Additionally, we observed an effective correlation between the small RNA-seq and q-PCR expression data for randomly selected miRNAs across five lactation stages. These results shown a similar expression pattern (average Person’s r = 0.71), which suggests that the sequencing data were reliable (Additional file 5: Figure S3).

Fig. 2
figure 2

Exosomal miRNAs expression profiles during lactation period. a Hierarchical clustering analysis for the normalized expression of conserved miRNAs among seven exosomal miRNA libraries based on Euclidean distance. b The number of expressed miRNAs at the two major stages and three early stages. c Top 10 miRNAs with the highest expression levels in exosomal miRNA libraries. Box Plot with different color represents the percentage of each miRNA expression in each library. Four co-expressed miRNAs in all 5 lactation stages are connected by lines (cli-miR-148a-3p, a potential biomarker for the quality control of mammalian milk can be seen marked with red line)

As shown in Fig. 2c, abundant miRNAs concentrated on few miRNAs of PM in five lactation stages. The top ten expressed miRNAs over five lactation stages correspond to 25 kinds of miRNAs and account for more than 69.48% (by total counts) of all the conserved miRNAs. Amid them, four conserved miRNAs (cli-miR-21-5p, cli-miR-148a-3p, cli-miR-10a-5p and cli-miR-26a-5p) were shared by all five lactation stages in top ten positions. These observations suggested that these miRNAs could be the important components of PM. Interestingly, based on annotation of the Pathway Central database (SABiosciences, MD, USA), a substantial proportion of miRNAs (including 3 of top ten co-expressed miRNAs) were designated as immune-related miRNAs. Similar to previous reports on mammalian breast milk [22,23,24], out of 139 immune-related miRNAs that are most relevant to the function of T-cell and B-cell activation, inflammatory response and autoimmunity, 58 (19.27%) are present and enriched in PM exosomal miRNA libraries (P < 10− 37, χ2 - test) (Fig. 3a). For instance, miR-146b, a negative regulator of the innate immune response by targeting NF-κB signaling [27]. miR-27b destabilizes peroxisome proliferator-activated receptor γ (PPARγ) mRNA, associating with chronic inflammatory diseases [28]. miR-200a-3p is associated with Hodgkin lymphoma [29], and potentially in regulating pathogen recognition in miiuy croaker through TLR1 targeting, an essential component of TLRs in bacterial pathogen defense [30]. The miR-26a, in transformed avian lymphocyte lines, regulates the expression of interleukin-2, which is a significant factor in development, differentiation, proliferation, and homeostasis of T cells [31]. Furthermore, we found that 38 (12.62%) growth-related miRNAs (P < 10− 29, χ2 - test) within cell differentiation and development pathways were also present in PM (Fig. 3b). These included miR-21-5p, miR-26a-5p, miR-10a-5p. Amid these, miR-21 has been reported to increase adipogenic differentiation of human adipose tissue-derived mesenchymal stem cells (hASCs) through modulation of TGF-β signaling pathway [32]. miR-26a regulates tissue and cell growth and differentiation [33, 34]. miR-10a belongs to the miR-10 family, a miRNA largely involved in development by targeting Hox transcripts in several species [35].

Fig. 3
figure 3

The distribution of immune and growth-related miRNAs in PM. a The accumulative perception of immune-related miRNAs in five lactation stages. b The accumulative perception of growth-related miRNAs five lactation stages. χ2-test: Number of immune and growth-related miRNAs and others detected in PM exosomes compared with the total entries in miRBase 21.0. *P < 0.05, **P < 0.01 vs human immune and growth-related miRNAs based on annotation in the Pathway Central database (SABiosciences, MD, USA)

Comparative miRNomes between PM and mammalian milk

To further enlighten the miRNAs in PM, we identified and compared the conserved orthologous miRNAs in PM and mammalian milk (Additional file 4: Table S2). As depicted in Fig. 4a, out of 301 mature miRNAs in PM, 41 orthologous miRNAs group (giving rise to 81 mature miRNA) were identified across all four species (Additional file 6: Table S3). To examine expression differences among the species, 41 orthologous miRNAs group were further analyzed. The results demonstrated that the libraries of porcine and bovine breast milk were clustered together, and then converged with human breast milk, whilst the libraries of PM constructed a completely separate branch (Fig. 4b). Unsurprisingly, the libraries clustered together based on their respective species of origin. PCA exhibited the specific regions for each species, and the variations in the conserved miRNA expression differed primarily by species (Fig. 4c).

Fig. 4
figure 4

Comparative analysis of conserved miRNAs in PM, human, bovine and porcine breast milks. a The ortholog groups of miRNA. b Hierarchical clustering analysis of conserved miRNAs across 23 libraries. c Principal component analysis (PCA) of conserved miRNAs across all 23 libraries. d The pairwise correlation of conserved miRNAs expression among different species

The pairwise correlation of conserved miRNAs expression across the species is shown in Fig. 4d. The expression of conserved miRNAs between human and bovine breast milk exhibited the highest correlation (average Pearson’s r = 0.75) compared to each of other libraries, followed by human vs. porcine breast milk, and bovine vs. porcine breast milk. However, the expression correlation of conserved miRNAs between human, bovine, porcine breast milk and PM were 0.60, 0.62 and 0.50, respectively.

The distinct miRNA expression patterns of PM during different lactation stages

Exosomal miRNAs expression pattern across the lactation period in PM were assessed by the analysis of sequencing data using the Short Time-series Expression Miner software (STEM) algorithm. The STEM clustering tool assigned each miRNA to the model profile that most closely matched its temporal expression profile [36]. Three significant model profiles (profiles 1, 2 and 3) marked with colors from the 50 distinct expression patterns were generated (Fig. 5a and Additional file 7: Table S4). Profiles 1, 2 and 3 comprised 20, 55 and 12 miRNAs, respectively. Shown with the same color, profiles 1 and 2 were characterized to form a cluster of similar profiles. Interestingly, the expression level of all three profiles gradually increased from 1 to 5 d, peaking on day 5 and then altered. However, the expression level for profiles 1 and 2 subsequently decreased to 20 d and 15 d, respectively. While, the expression level of the profile 3 remained high till 10 d and then decreased to 15 d; then it increased again to 20 d. These results indicate that the expressed miRNAs within these profiles may play different roles for the development of squabs. Therefore, we carried out the prediction of potential target genes for these miRNAs. As depicted in Fig. 5b, the target genes of miRNAs in profile 1 and 2 are mainly involved in cell migration (GO: 0016477), liver development (GO: 0001889), skeletal muscle tissue development (GO: 0007519) and mitogen activated protein kinase (MAPK) signaling pathway (gga04010), Wnt signaling pathway (gga04310), insulin signaling pathway (gga04910). Other functional terms including focal adhesion (gga04510), regulation of actin cytoskeleton (gga04810), endocytosis (gga04144) were also identified. Meanwhile, the target genes of miRNAs in profile 3 were enriched in Golgi apparatus (GO: 0005794), neuron projection (GO: 0043005) and heart development (GO: 0007507). The above results further suggest that the exosomal miRNAs in PM carry out array of important function in development of squabs.

Fig. 5
figure 5

The miRNA expression analysis by STEM clustering. a The profiles that had a statistically significant P-value and a number of miRNAs assigned using the ‘log normalize data’ option, with all other settings set to the defaults. Statistically significant profiles that were similar formed a cluster of profiles, and were shown in different color(s). b Gene Ontology (GO) and KEGG pathway analysis of potentially targeted miRNAs

Discussion

To date, no studies have investigated miRNA expression profiles in PM. Hence, to the best of our knowledge, this is the first report to identify exosomal miRNAs in PM. Because PM is regarded as indispensable for the growth and development of squabs, numerous studies have explored its nutritional composition [8,9,10]. Although the physiological synthesis and secretion of PM and mammalian milk differ, both provide similar nourishment for the growth and development of offspring.

Previous studies have demonstrated the existence of many exosomal miRNAs in mammalian milk [22,23,24]. Interestingly, we observed that PM also contains large amounts of exosomes, similar to those of mammalian milk, in the present study. Furthermore, a considerable proportion of miRNAs found in mammalian milk were also present in PM, which is consistent with previous reports on breast milk from humans [22, 37], pigs [23, 38], cows [39], and giant pandas [24].

Of these, miR-148a-3p, a potential biomarker for the quality control of mammalian milk, showed relatively high expression levels (in the top 10) in PM across five lactation stages. This indicated that it may be stably expressed and evolutionarily conserved from mammalian milk. Similarly, other miRNAs such as miR-21 and miR-26a were also enriched in PM; these are regarded as potential biomarkers for the quality control of raw milk and other milk-related products [39]. Furthermore, expression of the conserved miRNAs between PM and breast milk from humans and bovine was positively correlated (Pearson’s r = 0.60 and 0.62, respectively), while the average correlation between mammalian breast milk groups was 0.68. Based on these findings, it is tempting to speculate that pigeon ‘lactation’ and mammalian lactation generate similar functional products although they evolved independently [1].

The process of PM production begins when the germinal cell layer of the crop rapidly proliferates, followed by desquamation of lipid-rich cornified epithelium and subsequent accumulation in the crop lumen [6]. Our results suggest that some exosomal miRNAs in PM have a regulatory role in the proliferation and desquamation of crop cells. For instance, miR-21 in PM maintained high expression levels during lactation periods. Previous studies showed that miR-21 is one of several oncogenic miRNAs that control tumor cell proliferation, migration, and invasion [40, 41]. miR-21 was also found to be up-regulated by replicative and stress-induced senescence, despite being described as oncogenic [42], which could be associated with the hypoxic stress caused by a lack of blood supply to the rapidly proliferating germinal cell layer of the pigeon crop [7].

miR-205a is another miRNA that was abundantly expressed in PM. Previous reports revealed that it promoted the migration of keratinocytes via lipid phosphatase SHIP2, and that it was involved in regulation of the actin cytoskeleton and cell migration [43]. Of note, the accumulation of neutral lipids in keratinocytes is a unique trait of avian species [44], and pigeons were shown to utilize this ability to produce lipid-rich PM in their crop [7]. Furthermore, miR-21, miR-26a, and miR-27b, which were enriched in PM in this study, were previously shown to be associated with the adipogenic differentiation of adipose cells [45, 46]. Moreover, miR-22-3p, miR-148a-3p, and miR-26a-5p were highly expressed in chicken hepatocytes and exhibited critical roles in lipid metabolism [47]. Taken together, these results suggest that the highly expressed miRNAs in PM are likely to be associated with various important physiological functions, and may play essential roles in the generation of PM.

Through the analysis of miRNA expression patterns in PM, we identified three significant model profiles. GO and KEGG analyses of the putative target genes of highly expressed miRNAs in these three profiles revealed that enriched pathways were associated with organ development, particularly that of the heart, liver, and nervous system, and growth-related pathways including insulin, MAPK, and Wnt signaling pathways. These findings suggest that exosomal miRNAs in PM promote the growth and development of squabs.

Previous studies on the composition of PM indicated that it comprises proteins and lipids [11]. Hu et al. reported that the insulin signaling pathway IRS1/Akt/TOR regulates the synthesis of crop ‘milk’ protein [48], which is consistent with our KEGG analysis results. Moreover, the MAPK and Wnt signaling pathways may act on the stem cell population that gives rise to the proliferative germinal epithelium, which subsequently differentiates into keratinocytes that eventually form PM [7]. Additionally, the fatty acid precursors of crop triglycerides are likely obtained through the blood supply from oxidized triglycerides, which are transported from the liver or adipose tissue as fatty acids or very-low density lipoproteins that enter the cell by endocytosis [7]. Our KEGG pathway analysis further revealed that the endocytosis pathway is also enriched in exosomal miRNAs of PM.

In summary, these results indicate that PM miRNAs with different expression patterns play important regulatory roles in the synthesis and secretion of PM and the subsequent development of squabs.

Conclusions

In the present study, we have comprehensively enlightened and elucidated the distribution and expression profiles of exosomal miRNAs in PM and identified conserved miRNAs across species by small RNA-seq analysis, expanding the scope of the resources available to PM. Notably, similar to mammalian milk, a significant proportion of immune and growth-related miRNAs were present and enriched in PM. Besides, we also identified certain highly expressed miRNAs that may be involved in biosynthesis of PM; it is tempting to infer that our findings, therefore, prelude further detailed investigations elucidating the underlying essential roles of these highly expressed miRNAs. Nevertheless, based on our results, it is fascinating to speculate that the exosomal miRNAs in PM may be another factor affecting the development and growth of squabs.

Methods

Pigeon ‘Milk’ (PM) collection

The healthy pigeon squabs (Columba livia) were obtained from a pigeon farm in Wenjiang, Chengdu, Sichuan, China. PM samples were collected from squabs at five time points (1, 5, 10, 15 and 20 days after birth). As soon as the squabs were fed with PM, they were immediately anaesthetized with ether to minimize distress and subsequently euthanized. In brief, the pigeon squabs were respectively transferred to inverted beakers with cotton soaked in ether for about 30s. Subsequently, the squabs were euthanized by cervical dislocation. After 75% alcohol disinfection, their crops were surgically slit open by an incision of 2 cm in length. PM samples were collected from the crop of squabs into 15 ml sterile tubes and stored at − 80 °C for further analyses.

Preparation of exosomes

For the preparation of exosomes, the serial centrifugation and ultracentrifugation procedures were conducted as previously described [38], with minor modifications. Briefly, the PM was purified and stirred in phosphate buffered saline (PBS). PM samples were cleared of large aggregates as well as cell debris by sequential centrifugation at 2000×g for 30 min at 4 °C. The supernatant was sequentially subjected to centrifugations at 12,000×g for 45 min at 4 °C to remove other detritus. Filtration procedures were applied with syringe filters (0.22 μm) before the supernatant was prepared by ultracentrifugation at 160,000×g for 90 min. The final fraction was resuspended in 250 μl of PBS.

Atomic force microscope (AFM)

To characterize the morphology of isolated exosomes, PMs were diluted 1: 100 in deionized water and absorbed to freshly cleaved mica sheets for 20 min. The surplus solution was removed through careful absorption with a filter paper, and the mica was further dried before detection. Surface morphology was examined under an atomic force microscope (Asylum Research MFP-3D-Bio, Digital Instruments Inc., Santa Barbara, CA) as described by Sharma et al. [49].

Western blotting

Pigeon crop tissue and PM exosomes were lysed by protein lysis buffer (RIPA, Beyotime, Shanghai, China) supplemented with a protease inhibitor (Pierce, Rockford, IL, USA). The protein concentration was quantified by the BCA method. The lysates were boiled in SDS loading buffer (Beyotime), and were electrophoresed in 10% SDS-polyacrylamide gels and then transferred onto PVDF membranes (Millipore, USA). After blocking in 5% skim milk, the membranes were incubated with rabbit anti-alpha Tubulin (1: 1000, Abcam) or mouse anti-CD63 (1: 1000, Abcam) overnight at 4 °C. Then the PVDF membranes were washed with Tris-buffered saline with Tween 20 (TBST) and hybridized with Horseradish peroxidase-conjugated anti-rabbit or anti-mouse IgG, which was used as a secondary antibody (diluted 1: 5000 in TBST) for 2 h at RT. After washing three times, the protein bands were visualized with chemiluminescence reagents (Bio-Bad Laboratories) and scanned using ChemiDOC XRS instrument (Bio-Rad Laboratories).

Small RNA library preparation and sequencing

Exosomes were isolated from PM by serial differential centrifugation as previously described and collected in 250 μl PBS. Total RNA was extracted using Trizol-LS (Invitrogen, CA, USA) as per manufacturer’s instructions. The quality of RNA was examined by 1% agarose gel electrophoresis and further determined by the Agilent 2100 Bioanalyzer with RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA). For each PM library, small RNA ranging from 15 to 36 nt in length was purified by polyacrylamide gel electrophoresis and ligated using proprietary adaptors. The modified small RNA was then reverse-transcribed to cDNA and amplified by PCR. Finally, the libraries were sequenced on an Illumina HiSeq 2500 platform. We also downloaded 16 breast milk libraries of three mammals from NCBI’s Gene Expression Omnibus (GEO) Fdatabase, including human (GSE32253, n = 4) [22], bovine (GSE55144, n = 4) [50] and pig (GSE36590, n = 8) [23].

Identification of pigeon miRNAs and STEM analysis

To identify the pigeon miRNAs, initial sequence was subjected to a series of stringent filters (such as removing low quality-reads, repeated sequences and adaptor sequences) and the output was called clean data. Filtered sequences were then mapped to pigeon reference genome (Columba livia, ColLiv2) (https://www.ncbi.nlm.nih.gov/assembly/GCA_001887795.1) with stringent criteria (0 mismatch in the first 18 bp) using Bowtie software [51]. Since no exact annotation of pigeon miRNAs was recorded in miRBase 21.0, we performed alignment between the extended genomic sequence of all mappable reads and mature miRNA sequences of chicken (Gallus gallus), zebra finch (Taeniopygia guttata) and other mammals allowing no mismatch. Likewise, novel miRNAs were further predicted using the miRDeep2 core algorithm. For each sample, counts were first normalized by the total count of mappable reads which is denoted as reads per million (RPM), then standardized by the formula log2 (RPM + 1), which allowed unbiased comparison among samples. For each miRNA, the normalized number of counts was compared between groups.

A time-series analysis of the miRNA expression data was performed using STEM software (Short Time-series Expression Miner) http://www.cs.cmu.edu/~jernst/stem/ [52]. By STEM analysis, each miRNA was assigned to the model profile most closely matched to its time series based on the correlation coefficient. STEM was run using the log normalize data option, with all other settings set to the defaults.

Identification of orthologous miRNAs

As miRNAs that share highly similar sequences are likely to have similar functions, orthology groups of multi-species miRNAs were mainly constructed based on their similarity over their entire length. The identification of orthologous procedures were conducted with reference literatures [53, 54]. At first, we performed an all-against-all Blast search of miRNA precursors detected in all samples consisting of PM and breast milk libraries derived from human, bovine and pig. Blast hits with > = 75% of both precursor aligned and an identity of > = 70% were retained. Secondly, the miRNAs paired through blast hits were then grouped into homologous families using hcluster_sg [55]. To further eliminate the duplication miRNA genes in each homologous family and get the 1–1-1 orthology, we brought in a genomic location-matching approach based on BED tools [56], in which only miRNAs with the maximal overlap of precursor genomic coordinates were kept for each orthologous family.

Prediction and functional annotation of miRNA target genes

The algorithm TargetScan (http://www.targetscan.org/vert_71/) was applied to annotate miRNA target genes [57]. Gene Ontology (GO) and KEGG pathway analyses of target genes were performed using DAVID bioinformatics resources (https://david.ncifcrf.gov/) [58], with the Chicken (Gallus gallus) genome was selected as the background reference, which is the most closely related species for pigeon of the currently available entries in DAVID database. DAVID gene ontology enrichment was regarded as significant if the Benjamini P-value was smaller than 0.05.

Q-PCR

The expression changes of randomly selected miRNAs were determined by q-PCR approach. The cDNA was synthesized from miRNAs using the Mir-X™ miRNA First Strand Synthesis Kit (Clontech, Dalian, China) according to the supplier’s protocol. Quantitative PCR was carried out using SYBR Premix Ex Taq II (Takara, Dalian, China) on the CFX96™ Real-Time PCR Detection System (Bio-Rad, CA, USA). U6 snRNA were simultaneously used as internal control genes. The miRNA forward primers were obtained commercially from BGI (Shenzhen, China) and the universal reverse primer for miRNAs was offered by the Mir-X™ miRNA First Strand Synthesis Kit (Clontech, Dalian, China). The primer sequences are provided in Additional file 8: Table S5. Each miRNA sample was analyzed in triplicate and the 2−ΔΔCt method was used to calculate the relative expression levels of objective miRNAs.

Statistical analysis

All data from q-PCR are presented as means ± SD. Pearson’s correlation was used to determine the relationship between the small RNA-seq and q-PCR approaches. The χ2-test was used to determine statistical significance. A value of P < 0.05 was considered as significant (* P < 0.05, ** P < 0.01).

Abbreviations

AFM:

Atomic force microscope

GALT:

Gut associated lymphoid tissue

GO:

Gene ontology

hASCs:

Human adipose tissue-derived mesenchymal stem cells

KEGG:

Kyoto encyclopedia of genes and genomes

PBS:

Phosphate buffered saline

PM:

Pigeon ‘Milk’

RPM:

Reads per million

STEM:

Short Time-series Expression Miner software

TBST:

Tris-buffered saline with tween 20

References

  1. Gillespie MJ, Stanley D, Chen H, Donald JA, Nicholas KR, Moore RJ, et al. Functional similarities between pigeon ‘Milk’ and mammalian Milk: induction of immune gene expression and modification of the microbiota. PLoS One. 2012;7(10):e48363.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Horseman ND, Buntin JD. Regulation of pigeon cropmilk secretion and parental behaviors by prolactin. Annu Rev Nutr. 1995;15(1):213–38.

    Article  CAS  PubMed  Google Scholar 

  3. Dumont JN. Prolactin-induced cytologic changes in the mucosa of the pigeon crop during crop-“milk” formation. Z Zellforsch Mikrosk Anat. 1965;68(6):755–82.

    Article  CAS  PubMed  Google Scholar 

  4. Weber W. Zur Histologie und Cytologie der Kropfmilchbildung der Taube. Z Zellforsch Mikrosk Anat. 1962;56(2):247–76.

    Article  Google Scholar 

  5. Litwer G. Die histologischen Veränderungen der Kropfwandung bei Tauben, zur Zeit der Bebrütung und Ausfütterung ihrer Jungen. Z Zellforsch Mikrosk Anat. 1926;3(4):695–722.

    Article  Google Scholar 

  6. Gillespie MJ, Crowley TM, Haring VR, Wilson SL, Harper JA, Payne JS, et al. Transcriptome analysis of pigeon milk production - role of cornification and triglyceride synthesis genes. BMC Genomics. 2013;14(1):169.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Gillespie MJ, Haring VR, Mccoll KA, Monaghan P, Donald JA, Nicholas KR, et al. Histological and global gene expression analysis of the 'lactating' pigeon crop. BMC Genomics. 2011;12(1):452.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Davies WL. The composition of the crop milk of pigeons. Biochem J. 1939;33(6):898–901.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Goudswaard J, van der Donk JA, van der Gaag I, Noordzij A. Peculiar IgA transfer in the pigeon from mother to squab. Dev Comp Immunol. 1979;3(2):307–19.

    Article  CAS  PubMed  Google Scholar 

  10. Shetty S, Salimath PV, Hegde SN. Carbohydrates of pigeon milk and their changes in the first week of secretion. Arch Int Physiol Biochim Biophys. 1994;102(5):277–80.

    CAS  PubMed  Google Scholar 

  11. Shetty S, Bharathi L, Shenoy KB, Hegde SN. Biochemical properties of pigeon milk and its effect on growth. J Comp Physiol B. 1992;162(7):632–6.

    Article  CAS  Google Scholar 

  12. Shetty S, Hegde SN, Bharathi L. Purification of a growth factor from pigeon milk. Bba-Gen Subjects. 1992;1117(2):193–8.

    Article  CAS  Google Scholar 

  13. Hegde SN. Composition of pigeon milk and its effect on growth in chicks. Indian J Exp Biol. 1973;11(3):238–9.

    CAS  PubMed  Google Scholar 

  14. Guareschi C. Necessity of maternal alimentary factors for the growth of young pigeons. Boll Soc Ital Biol Sper. 1936;11:411–2.

    CAS  Google Scholar 

  15. Kim VN, Han J, Siomi MC. Biogenesis of small RNAs in animals. Nat Rev Mol Cell Biol. 2009;10(2):126–39.

    Article  CAS  PubMed  Google Scholar 

  16. Suzuki HI, Miyazono K. Emerging complexity of microRNA generation cascades. J Biochem. 2011;149(1):15–25.

    Article  CAS  PubMed  Google Scholar 

  17. Xiao C, Rajewsky K. MicroRNA control in the immune system: basic principles. Cell. 2009;136(1):26–36.

    Article  CAS  PubMed  Google Scholar 

  18. Kosaka N, Izumi H, Sekine K, Ochiya T. microRNA as a new immune-regulatory agent in breast milk. Silence. 2010;1(1):7–7.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Lässer C, Alikhani VS, Ekström K, Eldh M, Paredes PT, Bossios A, et al. Human saliva, plasma and breast milk exosomes contain RNA: uptake by macrophages. J Transl Med. 2010;9(1):1–8.

    Google Scholar 

  20. Merchant ML, Powell DW, Wilkey DW, Cummins TD, Deegens JK, Rood IM, et al. Microfiltration isolation of human urinary exosomes for characterization by MS. Proteom Clin Appl. 2010;4(1):84–96.

    Article  CAS  Google Scholar 

  21. Tan A, Rajadas J, Seifalian AM. Exosomes as nano-theranostic delivery platforms for gene therapy. Adv Drug Deliver Rev. 2013;65(3):357–67.

    Article  CAS  Google Scholar 

  22. Zhou Q, Li M, Wang X, Li Q, Wang T, Zhu Q, et al. Immune-related microRNAs are abundant in breast milk exosomes. Int J Biol Sci. 2012;8(1):118–223.

    Article  CAS  PubMed  Google Scholar 

  23. Gu Y, Li M, Wang T, Liang Y, Zhong Z, Wang X, et al. Lactation-related microRNA expression profiles of porcine breast milk exosomes. PLoS One. 2012;7(8):e43691.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Ma J, Wang C, Long K, Zhang H, Zhang J, Jin L, et al. Exosomal microRNAs in giant panda (Ailuropoda melanoleuca) breast milk: potential maternal regulators for the development of newborn cubs. Sci Rep. 2017;7(1):3507.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Munagala R, Aqil F, Jeyabalan J, Gupta RC. Bovine milk-derived exosomes for drug delivery. Cancer Lett. 2016;371(1):48–61.

    Article  CAS  PubMed  Google Scholar 

  26. Huang Z, Jebb D, Teeling EC. Blood miRNomes and transcriptomes reveal novel longevity mechanisms in the long-lived bat, Myotis myotis. BMC Genomics. 2016;17(1):906.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Taganov KD, Boldin MP, Chang KJ, Baltimore D. NF-kappaB-dependent induction of microRNA miR-146, an inhibitor targeted to signaling proteins of innate immune responses. Proc Natl Acad Sci U S A. 2006;103(33):12481–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Jennewein C, von Knethen A, Schmid T, Brune B. MicroRNA-27b contributes to lipopolysaccharide-mediated peroxisome proliferator-activated receptor gamma (PPARgamma) mRNA destabilization. J Biol Chem. 2010;285(16):11846–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Navarro A, Gaya A, Martinez A, Urbano-Ispizua A, Pons A, Balagué O, et al. MicroRNA expression profiling in classic Hodgkin lymphoma. Blood. 2008;111(5):2825–32.

    Article  CAS  PubMed  Google Scholar 

  30. Wang Y, Xu G, Han J, Xu T. miR-200a-3p regulates TLR1 expression in bacterial challenged miiuy croaker. Dev Comp Immunol. 2016;63:181–6.

    Article  CAS  PubMed  Google Scholar 

  31. Xu H, Yao Y, Smith LP, Nair V. MicroRNA-26a-mediated regulation of interleukin-2 expression in transformed avian lymphocyte lines. Cancer Cell Int. 2010;10(1):15.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Kim YJ, Hwang SJ, Bae YC, Jung JS. MiR-21 regulates adipogenic differentiation through the modulation of TGF-β signaling in mesenchymal stem cells derived from human adipose tissue. Stem Cells. 2009;27(12):3093–102.

    CAS  PubMed  Google Scholar 

  33. Wong CF, Tellam RL. MicroRNA-26a targets the histone methyltransferase enhancer of Zeste homolog 2 during myogenesis. J Biol Chem. 2008;283(15):9836–43.

    Article  CAS  PubMed  Google Scholar 

  34. Luzi E, Marini F, Sala SC, Tognarini I, Galli G, Brandi ML. Osteogenic differentiation of human adipose tissue-derived stem cells is modulated by the miR-26a targeting of the SMAD1 transcription factor. J Bone Miner Res. 2008;23(2):287–95.

    Article  CAS  PubMed  Google Scholar 

  35. Lund AH. miR-10 in development and cancer. Cell Death Differ. 2010;17(2):209–14.

    Article  CAS  PubMed  Google Scholar 

  36. Ahmad CM, Omaruddin RA, Brumbaugh CD, Tariq MA, Nader P. Identification of radiation-induced microRNA transcriptome by next-generation massively parallel sequencing. J Radiat Res. 2013;54(5):808–22.

    Article  Google Scholar 

  37. Simpson MR, Brede G, Johansen J, Johnsen R, Storro O, Saetrom P, et al. Human breast Milk miRNA, maternal probiotic supplementation and atopic dermatitis in offspring. PLoS One. 2015;10(12):e0143496.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Chen T, Xi QY, Ye RS, Cheng X, Qi QE, Wang SB, et al. Exploration of microRNAs in porcine milk exosomes. BMC Genomics. 2014;15(1):100.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Chen X, Gao C, Li H, Huang L, Sun Q, Dong Y, et al. Identification and characterization of microRNAs in raw milk during different periods of lactation, commercial fluid, and powdered milk products. Cell Res. 2010;20(10):1128–37.

    Article  CAS  PubMed  Google Scholar 

  40. Yan LX, Huang XF, Shao Q, Huang MY, Deng L, Wu QL, et al. MicroRNA miR-21 overexpression in human breast cancer is associated with advanced clinical stage, lymph node metastasis and patient poor prognosis. RNA. 2008;14(11):2348–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kumarswamy R, Volkmann I, Thum T. Regulation and function of miRNA-21 in health and disease. RNA Biol. 2011;8(5):706–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Dellago H, Preschitz-Kammerhofer B, Terlecki-Zaniewicz L, Schreiner C, Fortschegger K, Chang MWF, et al. High levels of oncomiR-21 contribute to the senescence-induced growth arrest in normal human cells and its knock-down increases the replicative lifespan. Aging Cell. 2013;12(3):446–58.

    Article  CAS  PubMed  Google Scholar 

  43. Yu J, Peng H, Ruan Q, Fatima A, Getsios S, Lavker RM. MicroRNA-205 promotes keratinocyte migration via the lipid phosphatase SHIP2. FASEB J. 2010;24(10):3950–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Vanhoutteghem A, Londero T, Ghinea N, Djian P. Serial cultivation of chicken keratinocytes, a composite cell type that accumulates lipids and synthesizes a novel beta-keratin. Differentiation. 2004;72(4):123–37.

    Article  CAS  PubMed  Google Scholar 

  45. Yao J, Wang Y, Wang W, Wang N, Li H. Solexa sequencing analysis of chicken pre-adipocyte microRNAs. Biosci Biotechnol Biochem. 2014;75(1):54–61.

    Article  Google Scholar 

  46. Arias N, Aguirre L, Fernandez-Quintela A, Gonzalez M, Lasa A, Miranda J, et al. Erratum to: MicroRNAs involved in the browning process of adipocytes. J Physiol Biochem. 2016;72(3):523–4.

    Article  CAS  PubMed  Google Scholar 

  47. Hong L, Zheng M, Jia L, Li Y, Xu C, Wang T, et al. Systematic analysis of the regulatory functions of microRNAs in chicken hepatic lipid metabolism. Sci Rep. 2016;6:31766.

    Article  Google Scholar 

  48. Hu XC, Gao CQ, Wang XH, Yan HC, Chen ZS, Wang XQ. Crop milk protein is synthesised following activation of the IRS1/Akt/TOR signalling pathway in the domestic pigeon (Columba livia). Brit Poultry Sci. 2016;57(6):855–62.

    Article  CAS  Google Scholar 

  49. Sharma S, Rasool HI, Palanisamy V, Mathisen C, Schmidt M, Wong DT, et al. Structural-mechanical characterization of nanoparticle exosomes in human saliva, using correlative AFM, FESEM, and force spectroscopy. ACS Nano. 2010;4(4):1921–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Sun J, Aswath K, Schroeder SG, Lippolis JD, Reinhardt TA, Sonstegard TS. MicroRNA expression profiles of bovine milk exosomes in response to Staphylococcus aureus infection. BMC Genomics. 2015;16(1):806.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Langmead B: Aligning short sequencing reads with bowtie. Curr Protoc Bioinformatics. 2010. Chapter 11: Unit 11.7.

  52. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7(1):191.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Santpere G, Lopezvalenzuela M, Petitmarty N, Navarro A, Espinosaparrilla Y. Differences in molecular evolutionary rates among microRNAs in the human and chimpanzee genomes. BMC Genomics. 2016;17(1):528.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Meunier J, Lemoine F, Soumillon M, Liechti A, Weier M, Guschanski K, et al. Birth and expression evolution of mammalian microRNA genes. Genome Res. 2013;23(1):34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Ruan J, Li H, Chen Z, Coghlan A, Coin LJ, Guo Y, et al. TreeFam: 2008 update. Nucleic Acids Res. 2008;36(Suppl 1):D735–40.

    CAS  PubMed  Google Scholar 

  56. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are MicroRNA targets. Cell. 2005;120:15–20.

    Article  CAS  PubMed  Google Scholar 

  58. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank all contributors of the present study. We also thank Dr. Haifeng Liu and students Guilin Li, Huan Li for their extensive help in the experiment.

Funding

This work was supported by grants from the National Natural Science Foundation of China (31472081 and 31522055), the Application Basic Research Plan Project of Sichuan Province (2016JY0167), the Sichuan Province & Chinese Academy of Science & Technology Cooperation Project (2017JZ0025), the National Program for Support of Top-notch Young Professionals, the Young Scholars of the Yangtze River.

Availability of data and materials

The raw data we obtained from small RNA-seq have been deposited in NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/projects/geo/) and are accessible through GEO Super Series GSE107312. The transcript annotation information for the miRNAs have been provided in the (Additional file 6: Table S3).

Author information

Authors and Affiliations

Authors

Contributions

YM, XW, SF, IHQ, KL and ML conceived and designed the study and drafted the manuscript. YM, GL, CN, YL, JX, XL and YW performed the experiments. SF, SH, XW, DL, YH, QT, AJ, JM and LJ analyzed all the experiments. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Xun Wang or Mingzhou Li.

Ethics declarations

Ethics approval and consent to participate

Research procedures involving animals were performed according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and approved by the Institutional Animal Care and Use Committee in College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Agilent 2100 analysis of total RNA from exosomes in PM. (PDF 952 kb)

Additional file 2:

Table S1. The sequencing results. (XLSX 9 kb)

Additional file 3:

Figure S2. Length distribution and expression distribution of identified exosomal miRNAs in PM. (PDF 396 kb)

Additional file 4:

Table S2. The mature miRNA group identified in PM and other three mammals. (XLSX 645 kb)

Additional file 5:

Figure S3. Q-PCR validation of miRNAs. (PDF 388 kb)

Additional file 6:

Table S3. The conserved miRNA group identified in 4 species. (XLSX 155 kb)

Additional file 7:

Table S4. miRNAs belonging to three Profiles as assigned by the STEM clustering tool. (XLSX 20 kb)

Additional file 8:

Table S5. The Primer sequences of q-PCR experiments. (XLSX 9 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ma, Y., Feng, S., Wang, X. et al. Exploration of exosomal microRNA expression profiles in pigeon ‘Milk’ during the lactation period. BMC Genomics 19, 828 (2018). https://doi.org/10.1186/s12864-018-5201-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5201-0

Keywords