Exosomal microRNAs in giant panda (Ailuropoda melanoleuca) breast milk: potential maternal regulators for the development of newborn cubs

The physiological role of miRNAs is widely understood to include fine-tuning the post-transcriptional regulation of a wide array of biological processes. Extensive studies have indicated that exosomal miRNAs in the bodily fluids of various organisms can be transferred between living cells for the delivery of gene silencing signals. Here, we illustrated the expression characteristics of exosomal miRNAs in giant panda breast milk during distinct lactation periods and highlighted the enrichment of immune- and development-related endogenous miRNAs in colostral and mature giant panda milk. These miRNAs are stable, even under certain harsh conditions, via the protection of extracellular vesicles. These findings indicate that breast milk may facilitate the dietary intake of maternal miRNAs by infants for the regulation of postnatal development. We also detected exogenous plant miRNAs from the primary food source of the giant panda (bamboo) in the exosomes of giant panda breast milk that were associated with regulatory roles in basic metabolism and neuron development. This result suggested that dietary plant miRNAs are absorbed by host cells and subsequently secreted into bodily fluids as potential cross-kingdom regulators. In conclusion, exosomal miRNAs in giant panda breast milk may be crucial maternal regulators for the development of intrinsic ‘slink’ newborn cubs.


Results and Discussion
Identification of endogenous exosomal miRNAs in giant panda breast milk. Exosomes were isolated from giant panda breast milk by ultracentrifugation and were investigated using atomic force microscope (AFM) at the nanometer scale. The results are indicative of that giant panda breast milk contains a large number of round or oval membrane vesicles, which had an approximate width of 130 nm and a height of 15 nm (Fig. 1B); they exhibited classic morphological ultrastructure similar to the exosomes derived from human saliva, urine and in vitro cells [21][22][23][24][25] . According to the recommendation of the International Society for Extracellular Vesicles (ISEV) 26 , the identity of milk-derived exosomes were also confirmed by Western blot analysis, which showing expected selective enrichment of exosomal protein markers CD63 and TSG101, and the absence of Tubulin in exosomes compared to the cellular fraction lysate of the breast milk (Fig. 1C). We also observed that the milk exosomes contained a considerable number of RNA sequences shorter than 100 nt in length ( Fig. 1D and We then used small RNA-seq to reveal the miRNA species and abundance in seven small RNA libraries isolated from the exosomes of five lactation periods. Each small RNA library generated ~9.50 million (M) (9.50 ± 1.31 M, n = 7) raw reads, resulting in more than 66.47 M raw reads for seven libraries. On average, 65% (65.55 ± 11.94%, n = 7) of the raw reads in each library passed the quality filter and were defined as reliable miRNA candidates ( Supplementary Fig. S2A), consistent with the canonical size range of mammalian miR-NAs ( Supplementary Fig. S2B), confirming the presence of miRNAs in the exosomes of giant panda milk. After mapping to the giant panda genome (Ailuropoda melanoleuca, AilMel1.0) and the known miRNA precursors (pre-miRNAs) of 24 mammals deposited in miRBase 21.0, a total of 1,103 pre-miRNAs encoding 1,236 mature miRNAs were identified in the seven libraries, including 320 conserved and 916 novel miRNAs (lineage specific) (Fig. 1E and Supplementary Table S1). There were distinct pre-miRNAs encoding identical mature miRNAs, resulting in 1,236 mature miRNAs corresponding to 1,191 unique miRNA sequences (Supplementary Table S2). The highly correlated biological replicates had similar frequency distributions for global miRNA expression (mean Pearson's r = 0.84) ( Supplementary Fig. S2C), suggesting the experimental reliability of the small RNA-sequencing data obtained in the present study. In addition, the clone sequencing and qRT-PCR analysis of 14 randomly selected miRNAs further confirmed the small RNA-seq results ( Supplementary Fig. S3).

Development-and immune-related miRNAs were enriched in giant panda breast milk.
To investigate the potential function of endogenous exosomal miRNAs, their global expression patterns were analyzed via hierarchical clustering. The miRNA expression profiles showed an obvious lactation-specific pattern. Two major branches were defined: one representing 0, 3 and 7 days and one representing 15 and 30 days ( Fig. 2A). This clustering pattern corresponded to the intrinsic differences of biochemical (e.g. vitamin, sugar, sodium and creatinine) and bioactive components (e.g. Immunoglobulin, interleukin and chemokine) between the giant panda colostrum (0 to 7 days after delivery) and mature milk (8 days after delivery) 27 , implying similar regulatory roles of exosomal miRNAs as potential functional ingredients in giant panda milk. Notably, we observed that the majority of abundant miRNAs were from a few miRNA species. As shown in Fig. 2B, the top 10 unique miRNAs with the highest expression levels accounted for more than 57.29% (57.29% to 69.08%) of the total read counts of all miRNAs for each library. The unified set of the top 10 unique miRNAs over the five lactation stages corresponded to 14 unique miRNAs, including six members of the let-7 family (i.e., ame-let-7a-1/-2/-3-5p, ame-let-7b-5p, ame-let-7c-5p, ame-let-7f-1/2-5p, ame-let-7g-5p and ame-let-7i-5p). In addition, seven of these top-ranked unique miRNAs (ame-let-7b-5p, ame-miR-92a-1/-2-3p, ame-miR-148a-3p, ame-miR-30a-5p, ame-let-7a-1/-2/-3-5p, ame-miR-181a-1/-2-5p and ame-let-7g-5p) were shared by all five lactation stages in the top 10 positions. Notably, most of the top-ranked miRNAs enriched in giant panda breast milk are also abundant in human 18 , porcine 28 and bovine 29 breast milk and have been associated with developmental processes and immune responses based on the annotation of the Pathway Central database (SABiosciences, MD, USA). For example, members of the let-7 family that are highly expressed in various species, including mammals, flies, worms, and plants 30 , are well-known regulators of development, cellular basal metabolism, and innate immune responses 31 . miR-148a regulates the innate response and antigen presentation of Toll-like receptor (TLR)-triggered dendritic cells by targeting CaMKIIα 32 . MiR-92a, member of the miR-17-92 cluster, negatively regulates TLR-triggered inflammatory responses in macrophages by targeting MKK4 kinase 33 . MiR-30a promotes B cell hyperactivity in patients with systemic lupus erythematosus by direct interactions with Lyn 34 . MiR-30d enhances invasion and immunosuppression during metastasis by targeting GalNAc transferases 35 . MiR-200a-3p is associated with various immunological diseases, such as Hodgkin lymphoma 36 and oral cancers 37 . Therefore, the high abundance of the above miRNAs in the milk of panda can affect the physiology, and particularly the immune system, of the newborn cubs. To induce an effect on the immune system of the cubs, the miRNA need to be horizontally transferred (or intestinally absorbed).

Distinct miRNA expression patterns during different lactation stages.
To visually illustrate the expression patterns of the miRNAs during different development stages, a short time-series expression miner (STEM) analysis was performed. The 1,194 unique miRNAs were clustered into 20 distinct expression patterns, of which 501 miRNAs were significantly over-enriched in three model profiles ( Fig. 3A and Supplementary  Table S3). There was a significantly larger-than-expected number of genes assigned to the model profile in the permutation runs (P < 0.01). The expression levels of miRNAs within a single profile are strongly correlated, whereas miRNAs belonging to different profiles generally showed no significant co-expression pattern (Fig. 3B), further confirming the results of STEM. Profile 1 comprised 334 unique miRNAs, whose expression gradually increased from 0 to 3 d, peaking on day 3 and subsequently decreasing to 15 d. Then, it increased again to 30 d. Profiles 2 and 3 comprised 114 and 54 unique miRNAs, respectively. The expression levels of miRNAs in profile 2 persistently increased following delivery, while the expression levels of profile 3 miRNAs decreased in this timeframe. These three profiles may indicate the different roles of corresponding miRNAs; therefore, we performed target prediction and functional annotation analyses of the conserved miRNAs in these three profiles using DIANA-mirPath 38 (see Methods for details). As shown in Fig. 3C, the target genes of conserved miRNAs that were highly expressed in the giant panda colostrum (0, 3 and 7 d; profiles 1 and 3) were primarily involved in immune and biosynthetic processes, including the "TLR signaling pathway, " "immune system process, " "proteoglycans in cancer, " "innate immune response, " "viral carcinogenesis" and "cellular nitrogen compound metabolic process". Recent reports have indicated that secreted exosomal miRNAs were enriched in the breast milk of various mammals, such as humans 8 , cows 19 , pigs 39 and rats 40 , and function as potential gene silencing signals between mothers and infants. The enrichment of these immune-related miRNAs in the giant panda colostrum may also imply their vital regulatory function in the neonatal immune system similar to other colostral immunological components. The colostral immunological components, such as immunoglobulin, cytokines, and leukocytes, that are transferred into the peripheral blood of the neonate play an important role in the development of the neonatal immune system 41 .
In contrast, the target genes of highly expressed miRNAs in mature milk (15 and 30 d; profiles 2) were primarily related to biosynthetic process, cell differentiation and gene expression, which are actively involved in body growth and development, such as "cell cycle, " "generation of precursor of metabolites and energy, " "cellular nitrogen compound metabolic process, " and "protein complex", and may be essential for the rapid development and growth of the premature newborn cubs. For example, miR-450a-5p promotes adipogenic differentiation by targeting WISP2 42 . MiR-34a promotes neural progenitor cell differentiation by targeting Numblike. MiR-135b-5p 43 and miR-664b-3p 44 can promote cell proliferation by targeting MEF2C and PLP2, respectively. MiR-183-5p 45 and miR-200b-3p 46 can positively regulate the development of inner ear. From 140 g (newborn) to 40 kg (weaning) of body weight, the growth and development of the giant panda primarily relies on the postnatal breast-feeding milk until weaning at 1-1.5 years old 3 based on many energy, protein and regulatory factors. This dependence may be the primary reason that miRNAs that play vital roles in growth and development are highly expressed in mature milk.
Resistance of milk exosomal miRNAs to the simulated gastrointestinal environment. The resistance to harsh conditions of the gastrointestinal environment is a precondition of maternal exosomal miR-NAs absorbed by the infant intestinal wall. To investigate the stability of exosomal miRNAs of giant panda breast milk, three endogenous and two exogenous miRNAs were randomly selected and subjected to a simulated gastrointestinal environment. Two exogenous synthetic worm-specific miRNAs (cel-miR-2 and cel-miR-93) were sharply degraded; however, the endogenous milk-derived miRNAs exhibited resistance to degradation when the milk was subjected to prolonged incubation at 37 °C, low pH or even treatment with high concentrations of exogenous RNase (Fig. 4). These results suggest that milk-derived miRNAs exist in a remarkably stable form, protected from endogenous RNase activity across various harsh conditions. This result is consistent with previous studies on blood 47 and human 8 and bovine 19 breast milk. The remarkable stability of the endogenous milk-derived miRNAs provides a basic condition for the dietary intake of maternal miRNAs by infants.
Exogenous plant miRNAs were detected in the milk exosomes of the giant panda. Interestingly, in the present study, we observed that the breast milk of the giant panda was colored yellowish green in 0 day colostrum and faded to milk white in 3 days; then, it gradually turned back to light yellowish green at 7 days (Fig. 5A), indicating that there are dynamic plant ingredients in the milk. Considering that giant pandas live primarily on bamboo and recent findings that plant miRNAs exist in mammalian body fluids 48 , we proposed that food-derived plant miRNAs also exist in breast milk. To test this hypothesis, bamboo miRNA species and their abundance were detected by mapping the small RNA-seq data to the bamboo genome and known plant precursor miRNAs in the miRBase. In total, 37 mature bamboo miRNAs corresponding to 27 precursor miRNAs were observed in the giant panda breast milk (Supplementary Table S4). To determine whether the bamboo miRNAs identified in giant panda breast milk are genuine plant miRNAs, we treated the total small RNAs isolated from giant panda milk with sodium periodate (oxidizing agent) 49 . Consequently, compared with the complete degradation of endogenous giant panda miRNAs (ame-let-7b-5p, ame-miR-181 and ame-miR-451), the bamboo-derived miRNAs (dla-miR-172d-3p, dla-miR-535-5p, dla-miR-1310-3p, dla-miR-1561-5p and dla-miR-2916-5p) exhibited similar resistance as synthetic miRNAs with 2′-Omethylated 3′ ends, highlighting the identity of bona fide plant miRNAs (Fig. 5B and Supplementary Fig. S4). Notably, the abundance levels of exogenous bamboo miRNAs in giant panda milk were significantly correlated with their intrinsic expression levels in bamboo leaf (Spearman's r = 0.58, P = 1.86 × 10 −4 ) (Fig. 5C, Supplementary Table 5) detected using a small RNA-seq approach. The moderate rather than high correlation may result from the different biodegradation degrees of bamboo miRNAs in long gastrointestinal tract before being absorbed, and from the selective packaging of intracellular miRNAs into exosomes/microparticles [50][51][52][53][54][55][56] . In addition, both small RNA-seq and stem-loop qRT-PCR results indicated that the abundance of exogenous bamboo miRNAs changed with the increasing greenness of giant panda milk (Fig. 5D,E), which may be highly associated with the stress fasting of the giant panda within the first three days after delivery 57 . This result emphasizes the potential oral intake of exogenous plant miRNAs from food sources in the giant panda. Functional enrichment analyses of the target genes indicated that the top ten exogenous plant miRNAs were primarily involved in basic cell metabolism, such as "protein binding, " "transcriptional activator activity, " "ubiquitin protein transferase activity" and "chromatin binding, " and neurodevelopment, such as "synapse organization, " "neuron migration" and "axon guidance" (Fig. 5F and Supplementary Table S6). These results indicated the potential role of exosome-loaded exogenous bamboo miRNAs in the postnatal development of giant panda cubs.

Materials and Methods
Animal ethics statement. All

Preparation of experimental animals and samples.
Three female giant pandas were sampled over multiple lactations, including 0, 3, 7, 15 and 30 days after delivery. Breast milk samples (5-10 ml) were collected from each stage, and a portion of each milk sample was stored at 4 °C for atomic force microscope (AFM) analysis. The remaining portion of each milk sample was centrifuged at 2500 × g for 10 min, twice, to remove cells and fat globules. The supernatant was centrifuged at 12,000 × g for 30 min at 4 °C to remove any cellular debris. The supernatants were then frozen in liquid nitrogen and subsequently stored at −80 °C until exosome enrichment and RNA extraction.
Atomic force microscope (AFM) assay. The breast milk exosomes for AFM analysis were prepared as described by Gu et al. 28 , with minor modifications. Briefly, 5 ml of breast milk was centrifuged at 2500 × g for 10 min, twice, to remove cells and fat globules. The supernatant was centrifuged at 12,000 × g for 30 min at 4 °C to remove any cellular debris. The supernatants were subsequently centrifuged at 120,000 × g for 4 hours, and the pellet was re-suspended in 10 ml of phosphate buffered saline (PBS). The exosome solutions were diluted 1:500 in deionized water and adsorbed to freshly cleaved mica sheets for 10 min. The redundant solution was removed by careful absorption with a filter paper, and the mica was further dried prior to detection. Surface morphology was examined under an AFM (Asylum Research MFP-3D-Bio, Digital Instruments Inc., Santa Barbara, CA) using identical conditions. The size of the exosomes was determined from a line profile of the AFM image.
Enrichment of exosomal RNA from breast milk. RNA in the breast milk exosomes was isolated as previously described 28 . In brief, raw milk (1.5 ml) was centrifuged at 2,000 × g for 10 min to remove fat globules. The supernatant was centrifuged at 12,000 × g for 30 min and further filtered through a 0.45 mm PVDF filter to eliminate cells and cellular debris. Approximately 1 ml of the supernatant was mixed with 500 ml of ExoQuick Exosome Precipitation Solution (SBI, CA, USA) and incubated at 4 °C for 12 hours. The ExoQuick/supernatant mixture was centrifuged at 1,500 × g for 30 min to obtain a beige exosome pellet, which was then re-suspended in 250 ml nuclease free water. Total RNA from the exosomes was extracted using TRIzol-LS (Invitrogen, CA, USA), following the manufacturer's instructions. Small RNAs were analyzed with the Agilent Bioanalyzer 2100 and the RNA 6000 Nano LabChip Kit (Agilent, CA, USA). Western blot analysis. Thirty micrograms of protein was separated on an 8% SDSPAGE gel and transferred to polyvinylidenedifluoride membranes (Millipore). Membranes were blocked and then incubated for 2 hours with either rabbit anti-alpha Tubulin (1:1000, Abcam), mouse anti-CD63 (1:1000, Abcam) or mouse anti-TSG101 (1:200, Abcam). Horseradish peroxidase-conjugated anti-mouse or anti-rabbit IgG was used as a secondary antibody (diluted 1:5,000 in TBST). Bands were scanned using a densitometer (GS-700; Bio-Rad Laboratories).
Small RNA sequencing. Total RNAs isolated from three individuals at four lactation stages (0, 7, 15 and 30 days after delivery) were pooled at equal quantities for library construction of each stage. To evaluate the variation among individuals, the total RNA from the milk-derived exosmes of three giant pandas at 3 days were used separately for library construction. In total, seven libraries were subjected to single-end sequencing of 36 nt reads using an Illumina Genome Analyzer II. The bioinformatics pipeline for miRNA discovery and profiling was used as previously described, with some improvements 58 . All reads were counted, and the identical reads were combined into a single type. The raw reads were subjected to a series of additional strict filters (such as the digital filters of base-call quality, read length, and sequence comparison) with acceptance criteria derived from the statistics of mammalian and plant miRNAs in miRBase 21.0. The reads originated from known RNAs of the giant panda and bamboo (i.e., mRNA; rRNA, tRNA, snRNA, and snoRNA), and repetitive sequence elements were also filtered. The raw reads were passed through a series of filters and mapped to all known pre-miRNAs of mammals and plants registered in miRBase 21.0, and the Giant Panda (Ailuropoda melanoleuca, AilMel1.0) and bamboo (http://www.ncgr.ac.cn/bamboo) genomes using NCBI Local BLAST. All candidates who did not map to any known pre-miRNAs but map to genome were considered as potential novel miRNAs. To further identify these potential novel miRNA candidates, MIREAP software 59 was used to predict novel miRNA by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the annotated small RNAs.
Short time-series expression miner (STEM) analysis. STEM analysis was used to visualize expression patterns for unique miRNAs. Each miRNA was assigned to the model profile most closely matched to its time series based on the correlation coefficient. The number of miRNAs assigned to each model profile was subsequently computed. The number of miRNAs assigned to a profile was estimated by randomly permuting the original time point values, renormalizing the miRNA expression values, assigning the miRNAs to their most closely matching model profiles, and repeating this process for a large number of permutations. The average number of miRNAs assigned to a model profile over all permutations was used as the estimate of the expected number of miRNAs assigned to the profile. The significance of the number of miRNAs assigned to each profile compared with the expected number was also computed 60 .

Databases of known signaling pathways. The pathways knowledge base developed by SABiosciences
(SABiosciences, MD, USA) (http://www.sabiosciences.com/pathwaycentral.php) was used to annotate the miR-NAs of interest in known pathways, as described previously 61 . Prediction and functional annotation of miRNA target genes. DIANA-mirPath 38 , a web-based computational tool, was employed to identify the target genes of giant panda miRNAs and the molecular Gene Ontology (GO) and pathways potentially altered by the expression of multiple miRNAs (http://www.microrna.gr/ miRPathv3). TargetScan human 5.2 62 was used for the target prediction of exogenous bamboo miRNAs based on seed sequences. Subsequently, the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms enriched in the predicted target genes were determined using Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources 63 . The predictions were made according to the interactions of human mRNA-miRNA due to the absence of giant panda and bamboo miRNAs in the current version of the above-mentioned algorithm.
qRT-PCR and clone sequencing validation of conserved and novel miRNAs. The changes in abundance of seven conserved and eight novel giant panda miRNAs over five lactation stages were determined by an EvaGreen-based qRT-PCR approach using a High-Specificity miRNA qRT-PCR Detection Kit (Stratagene, La Jolla, USA) on the CFX96 TM Real-Time PCR Detection System (Bio-Rad, CA, USA). The U6 snRNA, 5 S rRNA, and Met-tRNA were simultaneously used as endogenous control genes. All measurements contained a negative control (no cDNA template), and each RNA sample was analyzed in triplicate. Relative expression levels of objective miRNAs were calculated using the 2 −ΔΔCt method. To ensure the accuracy of the 2 −ΔΔCt results, we have carried out the optimization of annealing temperature by thermal gradient and evaluation of amplification efficiency in CFX96 TM Real-Time PCR Detection System (Bio-Rad, CA, USA), and all wells used in the subsequent analysis were met 0.95 < AE < 1.05, including the target genes and endogenous controls.
The PCR products of novel miRNAs were then cloned into plasmid using the TOPO TA Cloning Kit (Invitrogen, Paisley, UK), and subsequently sequenced. Stem-loop qRT-PCR was used to evaluate the exogenous plant miRNA level by Bulge-LoopTM miRNA qRT-PCR Starter Kit (RIBOBIO, China) according to manufacturer's instructions, and the transcriptional and qPCR primers were also synthesized by RIBOBIO (Guagnzhou, China) and shown in Supplementary Table S7. All measurements contained a negative control (no cDNA template), and each RNA sample was analyzed in triplicate.
Simulated gastrointestinal harsh treatments of miRNAs. Two exogenous miRNAs (cel-miR-2-3p, and cel-miR-39-5p of C. elegans), which have no sequence similarity to the giant panda miRNAs, were chemically synthesized. To determine the resistance and stability of endogenous exosomal miRNAs in the giant panda milk, exogenous miRNAs were added directly to raw milk, which were then subjected to three types of treatments: (A) incubation at 3 7 °C for 1, 2, or 3 hours; (B) treatment in pH 2.0 solution for 5, 30, 60, or 120 minutes; (C) treatment with RNase A (0.16 mg mL −1 ) and RNase T1 (0.4 U mL −1 ) (Fermentas, Shenzhen, China) for 5, 30, 60 or Scientific RepoRts | 7: 3507 | DOI:10.1038/s41598-017-03707-8 120 minutes at 37 °C. The treated milk samples was centrifuged at 2500 × g for 10 min, twice, to remove cells and fat globules. The supernatant was centrifuged at 12,000 × g for 30 min at 4 °C to remove any cellular debris. The supernatants were then used for total RNA extraction using TRIzol-LS (Invitrogen, CA, USA).
Statistical analysis. All data from qRT-PCR are expressed as mean ± SD. The statistical analysis was performed using SPSS 19.0 software. Equality of group variances were assessed by Kolmogorov-Smirnov test. Statistical significance was calculated by one way analysis of variance (ANOVA) with Tukey's post-hoc test for multiple groups or Student's t-test for comparisons of two groups. P < 0.05 was considered statistically significant (*P < 0.05; **P < 0.01). Spearman's correlations were used to examine the relationship between miRNA expression levels detected by small RNA-sequencing and qRT-PCR, and between bamboo miRNA expression levels presented in bamboo leaf and milk-derived exosome.

Conclusions
The present study presents a comprehensive survey of the exosomal miRNA transcriptome in giant panda breast milk during five lactation stages, shedding new light on the gene regulation network during the development of newborn cubs. Immune-and development-related endogenous miRNAs were observed to be enriched in the exosomes of giant panda breast milk, which together with evidence that these exosomal miRNAs are resistant to harsh gastrointestinal conditions, suggested that these exosomal miRNAs may be genetic material transferable from mother to infant and are essential for the development of newborn cubs. In addition, these findings support the viewpoint that food-derived exogenous plant miRNAs can be taken up and presented in mammalian breast milk.