This article is mentioned in:

Introduction

Down syndrome (DS) is the most frequent human chromosomal disorder with a frequency of 1/400 conceptions and 1/1,000 births worldwide (1,2). Since the initial discovery of Lejeune et al (3), it is known that the presence of full or partial chromosome 21 (Hsa21) in three copies (trisomy 21) in the cells of the affected subjects is responsible for the typical features of DS, in particular intellectual disability (ID), cardiovascular defects (4,5) and craniofacial dysmorphism. Importantly, a highly restricted 'Down syndrome critical region' of 34 kb on distal 21q22.13 appears to be specifically duplicated in all individuals with DS (6,7). In addition, a recent study conducted on plasma and urine samples of individuals with DS demonstrated a systematic deviation of metabolites involved in central metabolic processes associated with mitochondrial metabolism including the Krebs cycle, glycolysis and oxidative phosphorylation in DS (8).

Despite intense efforts of the scientific community, the pathogenic mechanism linking chromosome 21 and DS is remains largely uncharacterized, and pharmacological therapies targeting genes located on Hsa21 have not been developed so far. Therefore, current research shifted to study the involvement of non-coding RNAs (e.g., microRNAs) in the process (7). microRNAs (miRNAs/miR) are small single stranded nucleic acids that regulate gene expression post-transcriptionally via binding to different mRNA targets, resulting in inhibition of mRNA translation (9,10). Increasing evidence indicates that miRNAs serve important roles in a large variety of biological processes including development, differentiation, proliferation and apoptosis (11). It is well documented that intracellularly produced miRNAs can be secreted in the extracellular milieu, mostly associated with Ago proteins, bound to lipoproteins or secreted in extracellular vesicles (EVs) (12-15). Due to their nature, the latter two may provide an enriched and preserved source of miRNAs (16-19), therefore they may be better associated with the disease state than total plasma analysis (20,21), allowing a deeper understanding of the ongoing pathological processes and complementing previous studies.

Intriguingly, several circulating miRNAs have been differentially retrieved from either placenta or plasma samples of euploid and subjects with DS (22,23). However, only a few of them (miR-99a, miR-125b, let-7c and miR-155) were mapped on Hsa21 (24). It should also be noted that in other studies no Hsa21-derived miRNAs were identified to be differentially expressed in association with DS (25,26). This variability may derive from the type of sample analyzed and the extraction procedure used to recover the nucleic acids.

The present study focused for the first time on a subset of miRNAs carried by circulating biogenic nanoparticles, namely EVs and high-density lipoproteins (HDLs). Biogenic nanoparticles formulations, also referred to as nanoparticle-enriched fraction (NEF), were obtained from the plasma of young subjects with DS and their miRNA expression profiles systematically compared with those of NEF from their non-trisomic siblings. Interestingly, three novel candidate miRNAs emerged from the present study (miR-16-5p, miR-144-3p and miR-99b-5p) that could help unraveling molecular aspects underlying the complex DS pathogenic phenotype.

Materials and methods

Subject enrolment

For all participants involved in the present study (approved by the competent Ethics Committee of Sant'Orsola-Malpighi Hospital in Bologna, Italy; approval no. 39/2013/U/Tess) written informed consent was obtained from the parents in the case of subjects with age below 18 years and from the subjects themselves if aged over 18 years, according to the approved protocol.

Patient enrolment was performed at the Unit of Neonatology of Sant'Orsola-Malpighi Hospital in Bologna, Italy from February 2014 to March 2016. All methods were performed in accordance with ethical principles for Medical Research involving human subjects of the Helsinki Declaration.

The exclusion criteria of the enrolment to the study were: Age <2 years old, distress at birth, strongly premature (EG <35 weeks) and severe central nervous system (CNS) disease at birth. The presence or absence of each clinical feature and comorbidity was assessed based on the agreed clinical judgment of at least two pediatricians with long-time experience in the follow-up of children with DS, following a complete, systematic visit of the subject and an interview with the parents.

Case selection

In this study, a dataset containing 30 participants: 15 young subjects with DS (DS; 7 females and 8 males, aged 4-20 years, mean 11.71±5.12 years) and the corresponding 15 healthy siblings (C; 4 females and 11 males, aged 2-25 years, mean 14.24±6.42 years) was used. For each subject with DS, information was collected on the age of the parent, birth weight, development stage in the first months of life (sitting, babbling, walking and sphincter control), the presence/absence of 37 typical Down syndrome traits (dichotomous variables) and the evaluation of the neurological impairment (mild, moderate and severe). In the absence of the availability of cognitive test quantitative results, the latter evaluation was based, following accurate visit of the subject and an interview with the parents, on the agreed clinical judgment of at least two pediatricians with long-time experience in the follow-up of children with DS, aimed to assess if each subject fell in the relative top, medium or lower class of cognitive and global performance within the population.

Collection of blood samples

Blood samples (3 ml) were collected in EDTA-coated blood collection tubes and processed within 2 h from blood draw (27,28). All traceable identifiers were removed prior to analysis to protect patient confidentiality, all samples were analyzed anonymously (29). The sample was transferred in a new tube and centrifuged at 1,250 × g for 10 min at room temperature (centrifuge ALC 4235 A, rotor ALC T111). The plasma fraction was isolated and centrifuged for a second time at 800 × g for 30 min at room temperature (centrifuge ALC4214, rotor 6642) and the supernatant was transferred to new tubes without touching the pellet or the bottom of the tube, and divided in aliquots of 300 µl. All plasma samples were rapidly frozen and stored at −80°C until needed for subsequent analysis.

The exclusion criteria of plasma samples from the subsequent analysis were: Blood sample treatment after 2 h from the draw, or evident contamination of plasma samples by residual erythrocytes at the end of the treatment.

NEF separation from plasma and miRNA isolation

A total of 300 µl of plasma were used to obtain the NEF with Exosome Precipitation Solution (Serum/Plasma; Macherey-Nagel GmbH, Düren, Germany), according to manufacturer's protocol. Briefly the plasma was centrifuged at 10,000 × g for 15 min at 4°C to pellet intact cells, cellular debris and bigger vesicles. Then the supernatant was mixed with 'Exosome precipitation solution', incubated for 30 min on ice and finally centrifuged at 600 × g for 5 min at 4°C. The supernatant was kept for further analyses. The pellet was then resuspended in 300 µl RNAse free water and miRNAs were extracted using 'NucleoSpin miRNA Plasma' kit (Machery-Nagel GmbH), according to the manufacturer's protocol.

For reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analyses, 3 µl of a 4.16 nM solution of the synthetic miRNA cel-miR-39-3p from C. elegans (custom synthesized by Integrated DNA Technologies, Coralville, IO, USA) were added. miRNA fractions were eluted in 35 µl of nuclease-free water.

Western blotting

Proteins from NEF preparations obtained from C and DS plasma and associated supernatants were extracted with reducing SDS sample buffer (80 mM Tris, pH 6.8, 2% SDS, 7.5% glycerol, and 0.01% bromophenol blue) supplemented with 2% 2-mercaptoethanol for 5 min at 95°C and quantified by a Bradford's assay. A total of 50 µg of the preparations were loaded onto a 10% acrylamide/bis-acrylamide gel, electrophoretically separated by SDS-PAGE and transferred to polyvinylidene difluoride membranes (30,31). Membranes were then blocked at 37°C for 1 h in 5% non-fat milk, 10 mM Tris-HCl pH 7.5, 100 mM NaCl, 0.1% Tween-20 and probed with the following primary antibodies (all at 1:500 dilution) at 4°C for 16 h: Mouse anti-Lamp1 (BD Biosciences, San Jose, CA, USA; cat. no. 611042), mouse anti-cluster of differentiation (CD)63 (EMD Millipore, Billerica, MA, USA; clone RFAC4, cat. no. CBL553), rabbit anti-Adam10 (Acris OriGene Technologies GmbH, Germany, clone 3C6, cat. no. AP05830PU-N), mouse anti-TSG101 (Santa Cruz Biotechnology, Inc., Dallas, TX, USA; clone C-2, cat. no. sc-7964), mouse anti-Flotillin2 (BD Biosciences, San Jose, CA, USA; cat. no. 610383), mouse anti-CD81 (Santa Cruz Biotechnology, Inc.; clone B11, cat. no. sc-7637), mouse anti GM130 (BD Biosciences, San Jose, CA, USA; clone 35/GM130, cat. no. 610822) (32), mouse anti-ApoAI (Thermo Fischer Scientific, Inc.; clone 19H20L19, cat. no. 701239), mouse anti-Ago2 (Origene Technologies, Inc., Rockville, MD, USA; cat. no. TA352430) and incubated in the presence of specific horseradish-peroxidase conjugated immunoglobulin (Ig)G, Rabbit anti-mouse (cat. no. A90-117P) or Goat anti-rabbit (cat. no. A120-101P) (Bethyl Laboratories, Inc, USA, 1:3,000 dilution) at room temperature for 1 h (33,34). Immunoreactive bands were identified using the Luminata classico detection system (EMD Millipore). Images were acquired using a G:Box Chemi XT Imaging system and quantified using Gene Tools software version 4.01 (Syngene, Frederick, MD, USA) (35).

Atomic force microscopy (AFM) imaging

AFM imaging was conducted as described (36,37). Briefly, NEF preparations were diluted 1:200 with deionized water. A total of 5-10 µl of samples were then spotted onto freshly cleaved mica sheets (Grade V-1, thickness 0.15 mm, size 10×10 mm). All mica substrates were dried at room temperature and analyzed using a Nanosurf NaioAFM (Nanosurf AG, Liestal, Switzerland), equipped with Tap190AI-G tips (Budget Sensors; Innovative Solutions Bulgaria Ltd., Sofia, Bulgaria). Images were captured in tapping mode; the scan size ranged from 0.5-15 mm; the scan speed ranged from 0.6 to 1.5 sec × line.

miRNA expression by qPCR

Mature hsa-miR-16-5p, hsa-miR-99b-5p and hsa-miR-144-3p were amplified by a two-step Taq-Man RT-PCR analysis, using primers and probes obtained from Thermo Fisher Scientific, Inc., (hsa-miR-144-3p, cat. no. 002676; hsa-miR-16-5p, cat. no. 000391; hsa-miR-99b-5p, cat. no. 000436; cel-miR-39-3p, cat. no. 000200). cDNA was synthesized from 5 µl of the each RNA fraction in 15 µl reactions, using TaqMan MicroRNA Reverse Transcription kit (Thermo Fisher Scientific, Inc.). The reverse transcription reaction was performed by incubating the samples at 16°C for 30 min, followed by incubation at 42°C for 30 min and 85°C for 5 min. The RT-qPCR reaction (20 µl) contained 1.3 µl of reverse transcription product, 10 µl of Taq-Man 2X Universal PCR Master Mix and 1 µl of the appropriate TaqMan MicroRNA Assay (20X) specific for the miRNA targeted by the assay. The PCR mixtures were incubated at 95°C for 10 min, followed by 40 cycles at 95°C for 15 sec and 60°C for 60 sec. PCRs were performed in triplicate using the 7500 real time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The expression of miRs was based on the 2−∆∆Cq method (38), using cel-miR-39-3p as reference.

Bioinformatic analysis

For functional enrichment analysis of gene ontology terms and molecular signalling pathway, the following online software packages were used: Database for Annotation, Visualization and Integrated Discovery (DAVID) (39) and WebGestalt (40) and Kyoto Encyclopedia of Genes and Genomes analysis (KEGG; www.genome.jp/kegg/). The human disease database MalaCards (41) was used to search the diseases and the genes associated with DS. The miRNA targets were predicted by miRWALK (42). Considering that these procedures may have a certain false positive rate, the target genes of miRNA only considered the genes predicted by at least 7 procedures.

Statistical analysis

Descriptive statistics

The descriptive statistics for all quantitative variables available in the dataset (minimum, maximum, first, second and third quartile, mean ± standard deviation, number and corresponding percentage of missing values) and the frequency distributions of 37 typical DS pathologies (presence/absence of the pathology-in absolute value and %-) have been computed.

The investigation on the linear association between the quantitative variables collected on 15 individuals with DS has been performed by calculating the Pearson correlation coefficient ρ.

For the graphical representation of the significant associations (correlation tests with P<0.05) between variables under inspection, bubble charts were used that introduce in a simple scatterplot a third dimension, namely the diameter of the bubble. This dimension is proportional to the categorical variable called 'number of comorbidities'. In detail, for each individual with DS his/her number of comorbidities were summed, obtaining a quantitative score; then the first (Q1) and third (Q3) quartile were computed and the quantitative score were turned into the following three categories: i) First quartile (Q1=17.5), each subject with number of comorbidities <17.5 was classified as Low comorbidities; ii) third quartile (Q3=21), each subject with number of comorbidities >=17.5 and <21 was classified as Medium comorbidities; and iii) each subject with number of comorbidities >=21 was classified as High comorbidities.

Consequently, the scatterplots that highlight the significative association between quantitative variables report the dot circumference proportional to the number of assigned comorbidities: Small dots when subjects with DS are classified as 'Low comorbidities' (green); medium size dots for 'Medium comorbidities' (orange) and big dots when they are labeled as 'High comorbidities' (red).

For each of the qualitative variables reported (Table SI; considering also the neurological severity), the median values of the miRNAs expression levels in correspondence of subgroups of subjects defined by its categories were computed. Then, the non-parametric Wilcoxon signed-rank test (or Kruskal Wallis test when the subgroups are >2) was computed in order to compare the medians obtained (Table SII). These non-parametric tests do not require any assumptions on the shape distributions of the three miRNAs and permit analysis of the possible association between clinical variables and miRNAs expression. The unique significant associations were visualized by means of a boxplot for each category (Low, Medium and High comorbidities).

Heatmap generation

The data matrix with the investigated miRNAs was visualized by means of a heatmap. Starting from the standardized matrix (each miRNA has mean 0 and standard deviation 1), a heatmap was generated to visualize the miRNA expression (low values are represented with blue while high values with red). In this graph, similar values were placed near each other according to the clustering algorithm used in the analysis, thereby realizing one dendrogram appended on the left y-axis.

k-medoids cluster analysis

Starting from the data matrix containing only miRNA expression levels, we use the k-medoids cluster analysis generating k=2 groups. For evaluating the distance between subjects the Manhattan distance was used, which is based on absolute value distance. This choice should provide more robust results, whereas Euclidean distance would be influenced by unusual values in the data matrix.

Inside each group, subjects exhibit a high degree of similarity based on miRNA expression while subjects belonging to different groups are as dissimilar as possible. The k-medoids cluster analysis (known also with the acronym PAM-Partitioning Around Medoids) searched for k representative subjects (in the present study k=2) among all the subjects in the dataset. These representative subjects are called medoids (43) and the algorithm assigns each patients of the dataset to nearest medoid.

By means of silhouettes a graphical representation of cluster analysis was provided where the entire clustering was displayed by plotting the silhouette into a diagram highlighting the cluster quality. Denoting Cluster 1 with C1 and Cluster 2 with C2, silhouette for the subject i contained in C1 was:

From the ratio above, it is evident that each silhouette belongs to the interval -1≤ s(i) ≤+1. This means that when s(i) is close to 1 the within dissimilarities a(i) is smaller than the between dissimilarities b(i) and i is well classified. When s(i) is about 0, it is not clear where i should be assigned (cluster C1 or C2) and finally when s(i) is close to -1, a(i) is larger than b(i) and this means that i has been misclassified.

All statistical analyses were performed using R version 3.4.4.

Results

NEF separation, characterization and miRNA extraction

NEF samples were obtained with a commercial polymer-based precipitation (PBP) kit as described in the Materials and methods section. PBP is a simple, timesaving and cheap EV concentration method (44) characterized by low specificity. Formulations obtained by PBP consist in a heterogeneous mix of nanosized particles that overlap EVs in size and morphology, including HDLs.

NEF morphological properties were determined by AFM (45). The preparation was adsorbed on a freshly cleaved mica surface and imaged under ambient conditions. Samples were composed by spherical nano-objects with a size ranging from 10s to a few 100s of nm, which are visible despite the presence of background material, very like the polymeric matrix on which the precipitation kit is based (44,46) (Fig. 1A). The resulting fraction is mainly composed by small EVs and HDLs, as can be deduced by the enrichment of specific EV markers (47), including Lamp1, CD63, Adam 10, TSG101, Ago2 and CD81, as compared with the supernatant fractions by western blot analysis (Fig. 1B). On the other hand, NEF is devoid of both GM130s and Flotillin, excluding the presence of intracellular contaminants and implying that that large EVs are not efficiently enriched by this protocol. On the other hand Apo-A1, a biomarker characteristic of HDL lipoproteins (48), is detectable in supernatant fractions and in a lesser extent also in NEFs (Fig. 1B). Overall, the above results indicate that the final NEF is specifically enriched with nanosized RNA-carriers, namely small EVs and HDLs.

Figure 1 NEF characterization by AFM and western blotting. (A) Morphological analysis by AFM of C and DS NEF preparations. Images are representative of 3 independent experiments. (B) NEF obtained from C and DS plasma and the corresponding supernatants (SNC and SNDS) were immunoblotted for the indicated extracellular vesicles and high density lipoprotoeins markers (see text for details). Due to antibody-unspecific signals, the original image of western blotting membranes has been cropped. The uncropped, original image is available upon request. DS, Down syndrome; C, healthy control; NEF, nanoparticle-enriched fraction; AFM, atomic force microscopy; miR, microRNA; SNC, SNDS, C and DS plasma corresponding supernatants; M, marker; CD, cluster of differentiation.

NEF miRNA extraction and expression analysis

miRNAs were subsequently extracted from NEF using commercial columns. First, Agilent miR microarray analysis aimed to search for NEF miRNA content and was performed on samples derived from 4 couples of siblings: 4 subjects with DS and their 4 correspondent healthy siblings. This allowed identification of the most abundant miRNAs in the NEF recognized at least by 3 probes and expressed in all the 8 subjects included in the analysis (Table SIII). Consequently, the expression levels of the three most abundant miRNAs (hsa-miR-16-5p, hsa-miR-144-3p and hsa-miR-99b-5p) were further assessed in a larger cohort of DS and C by qPCR (15 couples of siblings; n=30).

As presented in Fig. 2A-C, the expression levels of these three miRNAs were very different among the couples of siblings considered and in general the results demonstrated an increased average expression trend in DS NEF with respect to C [miR-16-5p: Mean Relative Quantification (RQ) of DS: 29.5±13.3, mean RQ of C: 11.9±3.3, Fig. 2D; miR-144-3p: Mean RQ of DS 66.6±18.7, mean RQ of C: 49.5±14.3, Fig. 2E; miR-99b-5p: Mean RQ of DS: 42.0±15.8, mean RQ of C: 27.7±10.6 Fig. 2F].

Figure 2 Expression levels of miRNAs in NEF from the couples of siblings. Relative expression levels of (A) mature miR-16-5p, (B) miR-144-3p and (C) miR-99b-5p in each of the 15 siblings couples obtained by quantitative polymerase chain reaction. Average expression levels of (D) miR-16-5p, (E) miR-114-3p and (F) miR-99b-5p in individuals with DS and in control subjects. miR, microRNA; DS, Down syndrome; C, healthy control.

Since hsa-miR-16-5p, hsa-miR-144-3p and hsa-miR-99b-5p have all been previously identified as exosome-associated miRNAs and, in certain cases, associated with distinct pathologies (49-52), although never associated with the DS phenotype previously, it was decided to statistically correlate the expression levels of the three miRNAs with the available quantitative clinical data of the participants (quantitative variables, Table SIV) in order to have insight into their functional role in DS.

Correlation analysis of miRNA expression with clinical features of subjects with DS

In Table I the correlation matrix computed on the quantitative variables obtained from the subjects with DS was reported (correlation test with P<0.05). Bubble chart analysis (Fig. 3A-E) was subsequently performed only on the significant associations, introducing a third dimension, the diameter of the bubble. This dimension was proportional to the categorical variable named 'number of comorbidities', obtained as described in Material and methods section and reported in Table SV for each subject. An inverse significant linear association between the clinical feature 'development babbling' expressed in months and miR-16-5p and miR-144-3p expression levels was evident (Fig. 3A and B). From the bubble charts it is clear that these two miRNAs demonstrated high expression levels in subjects with DS with high number of comorbidities and beginning to babble earlier.

Figure 3 Bubble charts for the graphical representation of the significant association (ϱ > |0.7|, P<0.05) between the quantitative variables. The bubble diameter is proportional to the variable Number of comorbidities: High number of comorbidities is represented with a red circle while few numbers of comorbidities with a green circle and the orange circles are a middle way. (A) Linear negative association between miR-16-5p expression levels and the clinical feature 'development babbling' expressed in months. (B) Linear negative association between miR-144-3p and the clinical feature 'development babbling' expressed in months. (C) Linear positive association between miR16-5p and miR-144-3p expression levels. (D) Linear positive association between miR16-5p and miR-99b-5p expression levels. (E) Linear positive association between miR144-3p and miR-99b-5p expression levels.

Table I Correlation matrix.

Table I

Correlation matrix.

Variables	Mother age	Father age	Birth weight (g)	Sitting months	Babbling months	Walking months	Sphincter control months	miR-16-5p	miR-144-3p	miR-99b-5p
Mother age	1,0000
Father age	0,9197	1,0000
Birth weight (g)	−0,0982	−0,1577	1,0000
Sitting months	0,2095	0,2532	−0,4586	1,0000
Babbling months	0,1798	−0,0132	−0,1090	−0,0157	1,0000
Walking months	0,3488	0,1835	0,2421	0,2210	0,1501	1,0000
Sphincter control months	0,2018	0,2531	−0,3816	0,7177	0,2275	0,5274	1,0000
miR-16-5p	−0,1994	−0,0205	−0,0940	−0,1439	−0,7441	−0,4159	−0,2803	1,0000
miR-144-3p	−0,1070	0,0451	−0,2181	0,1122	−0,8070	−0,1911	−0,0212	0,8383	1,0000
miR-99b-5p	−0,1954	−0,0176	−0,2089	−0,1102	−0,6024	−0,1794	−0,0208	0,8780	0,8516	1,0000

[i] In bold correlation with P<0.05 (a test on the significance of the correlation coefficient was performed).

Concerning miR-99b-5p, it was observed that the median value of the expression level was significantly increased (Kruskal test: P=0.05) in individuals with DS presenting high comorbidities (Fig. 4 and Table SII) suggesting a potentially more relevant role of miR-99b-5p upregulation in this subgroup of subjects with DS. Instead, no significant differences were identified for miR-144-3p and miR-16-5p (Table SII).

Figure 4 Box plot on miR-99b-5p expression level using like grouping variable classification based on the number of comorbidities (Low, Medium and High). The individuals with Down syndrome bearing high comorbidities display the highest miR-99b-5p expression levels. miR, microRNA.

Fig. 3C highlights a significant linear positive association between miR16-5p and miR-144-3p expression. The same consideration can be extended to miR16-5p versus miR-99b-5p and miR144-3p versus miR-99b-5p, respectively (Fig. 3D and E). Furthermore, individuals with DS with the highest expression levels of miRNAs exhibited a high number of DS-associated typical clinical features. It is interesting to outline that when the same correlations were calculated by Pearson correlation analysis on the healthy siblings (on the basis of the descriptive statistics reported in Table SVI), we obtained low values of ρ (miR-16-5p vs. miR-144-3p: ρ= 0.1895; miR-16-5p vs. miR-99b-5p: ρ= 0.5123; miR-144-3p vs. miR-99b-5p: ρ=0.6820). This well agrees with the finding that the three selected miRNAs were upregulated in the DS cases in the present study, indicating a possible role in DS-associated typical clinical anomalies, probably mainly in subjects with 'high comorbidities'. As expected, by repeating the same analysis on all the 30 participants (DS+C), significant positive correlations were again obtained between miRNAs (Fig. 5A-C).

Figure 5 Scatter plots (with r and corresponding P-value) computed on all subjects (DS+C) evidencing the significant linear positive association between the miRs expression levels: (A) miR-16-5p vs. miR144-3p, (B) miR-16-5p vs. miR-99b-5p and (C) miR-144-3p vs. miR-99b-5p expression. Scatter plot (with ρ and corresponding P-value) of each miR expression levels [(D) miR-16-5p; (E) miR-144-3p and (F) miR-99b-5p], measured on DS vs. their corresponding siblings, C. The unique significant positive linear correlation is for miR-99b-5p. miR, microRNA; DS, Down syndrome; C, healthy.

ϱ was also computed on the expression level of each miRNA comparing couples of siblings (Fig. 5D-F). The only significant positive linear correlation was for miR-99b-5p (Fig. 5F). By computing ρ for couples of unrelated subjects (C subjects and subjects with DS were repeatedly randomized computing the correlation coefficients and corresponding test) it was noted that there was no correlation for miR-99b-5p (data and Figure provided upon request). This result demonstrates the same expression trend of miR-99b-5p among brothers, demonstrating the importance of the genetic background when analyzing miRs.

miRNA expression was visualized using heatmaps either on the entire sample (DS+C) (Fig. 6A), or separating C subjects and subjects with DS (Fig. 6B and C, respectively). The subjects with the highest miRNA expression (at the top three positions of the graph) belonged to the DS group. Furthermore, the C subjects at the top positions were siblings of the DS with highest miRNA expression. Focusing the attention on heatmaps in Fig. 6B and C, it was evident that miRNA expression levels were different between DS and C (in DS the color key is red-shifted compared with C); among the subjects with DS there was only one case exhibiting very low miRNA expression and notably this subject also presented low comorbidities; the miRNAs less expressed were miR-144-3p and miR-99b-5p in DS and C, respectively.

Figure 6 Heatmaps of the expression levels of miR-16-5p, miR-144-3p and miR-99b-5p. (A) Heatmaps on the entire sample DS+C, (B) on subjects with DS and (C) on healthy individuals. DS, Down syndrome; C, healthy individuals.

Using the Partitioning Around Medoids algorithm (PAM) 2 clusters of subjects were identified (Fig. 7A): The first one, called C1, contained subjects with low median and mean expression values for each miRNA with a prevalence of C subjects (55% of individuals in C1 were C which corresponded to 73% of healthy siblings in the study, see Table SVII for major details). Cluster 2, called C2, contained a number of subjects with DS with increased miRNA expression (60% of individuals in C2 were DS which corresponded to 40% of subjects with Down syndrome in the study). Furthermore, certain subjects with DS in C2 were classified together with their siblings (DS10 and C10, DS12 and C12) confirming similar miRNA expression levels.

Figure 7 Partitioning around medoids algorithm based on miRs expression levels of DS and C, highlighting 2 well separated clusters [Cluster 1(C1), subjects indicated with circles and cluster 2, (C2), subjects indicated with triangles]. (A) Silhouettes computed for each subject in the analysis. The ID of each subject is reported on the left side of the graph and the DS silhouettes are colored using the classification based on number of comorbidities. Cluster 1 (C1) contains subjects with low median and mean expression values for each miR with a prevalence of C subjects. Cluster 2 (C2) (comprising 60% of individuals with DS) with higher miR expressions and furthermore 60% of their corresponding siblings belongs to the same cluster. (B) miR, microRNA; DS, Down syndrome; C, healthy.

From the PAM algorithm it was also clear that miRNA levels were able to generate clusters of subjects that were well classified. In fact, the 'within dissimilarities' was smaller than the 'between dissimilarities' and 73% of individuals were well classified, as presented in Fig. 7B, where 22 subjects on 30 analyzed (on y-axis) exhibited a width silhouette s(i) (on x-axis) higher than 0.4. Therefore, this representation confirmed the goodness of the classification induced by the miRNAs matrix (note that s(i) reaches its maximum in correspondence of 1). Furthermore, in Fig. 7B each subject with DS is colored using the classification based on number of comorbidities: In C2 cluster, subjects with DS with high comorbidities prevailed.

miRNAs gene target prediction and enrichment study

The main putative targets of the 16 most abundant miRNAs identified in the NEFs (Table SIII) were subsequently predicted using the online software miRWalk (6). The analysis in the present study identified 2,515 targeted genes, including 25 of the 536 genes located on Hsa21. Functional annotation analysis of the 2,515 genes was performed next. In a GO analysis conducted by DAVID (1,865 unique IDs identified), the candidate genes were analyzed in all categories including biological process (Table II), molecular function (Table III) and cellular component (Table IV). For the biological process category, the first 4 enriched terms were referred to the transcriptional process; the fifth biological term was 'nervous system development' with 53 predicted genes included in this term (P=0.0140). In the molecular function category, 1,010 genes were included in 'protein binding' term and the other enriched terms belonged essentially to 'transcription process' categories. In the cellular component category, the 'nucleus' was the predominant term with 620 associated-genes, the fourth term was 'neuronal cell body' with 54 genes (P=0.0080). Therefore, >50 candidate miRNAs target genes may have a role in 'nervous system development' or in the 'neuronal cell body' formation.

Table II GO Top 10 significant terms belonging to the 'biological process' category.

Table II

GO Top 10 significant terms belonging to the 'biological process' category.

GO term (biological process)	Number of genes	%	P-value	Benjamini corrected P-value
Positive regulation of transcription. DNA-templated	102	5.5	1.40×10−11	6.80×10−08
Positive regulation of transcription from RNA polymerase II promoter	163	8.7	2.80×10−11	7.10×10−08
Negative regulation of transcription from RNA polymerase II promoter	126	6.8	2.50×10−10	4.20×10−07
Transcription from RNA polymerase II promoter	96	5.1	1.40×10−09	1.70×10−06
Nervous system development	53	2.8	1.40×10−05	1.40×10−02
Protein phosphorylation	73	3.9	4.90×10−05	4.00×10−02
Steroid hormone mediated signaling pathway	17	0.9	9.30×10−05	6.40×10−02
Dentate gyrus development	9	0.5	1.00×10−04	6.00×10−02
Stem cell differentiation	11	0.6	1.60×10−04	8.40×10−02
Transcription. DNA-templated	241	12.9	1.9×10−4	9.1×10−2

[i] Bold font indicates typical Down syndrome features. GO, Gene Ontology.

Table III GO Top 10 significant terms belonging to 'molecular function' category.

Table III

GO Top 10 significant terms belonging to 'molecular function' category.

GO term (molecular function)	Number of genes	%	P-value	Benjamini corrected P-value
Protein binding	1,010	54.2	1.6×10−13	2.4×10−10
Transcriptional activator activity. RNA polymerase II core promoter proximal region sequence-specific binding	56	3.0	1.0×10−9	7.5×10−7
Transcription factor activity. sequence-specific DNA binding	154	8.3	1.1×10−9	5.2×10−7
Sequence-specific DNA binding	88	4.7	5.1×10−7	1.9×10−4
RNA polymerase II core promoter proximal region	66	3.5	6.9×10−7	2.0×10−4
sequence-specific DNA binding
Protein serine/threonine kinase activity	66	3.5	5.2×10−6	1.3×10−3
RNA polymerase II regulatory region sequence-specific	41	2.2	2.9×10−5	6.0×10−3
DNA binding
GDP binding	17	0.9	4.2×10−5	7.7×10−3
Zinc ion binding	157	8.4	4.7×10−5	7.6×10−3
Transcriptional repressor activity. RNA polymerase II core promoter proximal region sequence-specific binding	26	1.4	6.6×10−5	9.5×10−3

[i] GO, Gene Ontology.

Table IV GO Top 10 significant terms belonging to 'cellular component' category.

Table IV

GO Top 10 significant terms belonging to 'cellular component' category.

GO term (cellular component)	Number of genes	%	P-value	Benjamini corrected P-value
Nucleus	620	33.2	1.1×10−7	8.2×10−5
Nucleoplasm	345	18.5	2.2×10−7	8.0×10−5
Transcription factor complex	38	2.0	3.8×10−5	9.2×10−3
Neuronal cell body	54	2.9	4.3×10−5	8.0×10−3
Cytosol	383	20.5	4.6×10−5	6.8×10−3
Protein-DNA complex	11	0.6	5.7×10−5	7.1×10−3
Perinuclear region of cytoplasm	89	4.8	1.5×10−4	1.6×10−2
Golgi apparatus	116	6.2	2.1×10−4	1.9×10−2
Golgi membrane	84	4.5	3.1×10−4	2.5×10−2
Cytoplasm	567	30.4	3.6×10−4	2.6×10−2

[i] Bold font indicates typical Down syndrome features. GO, Gene Ontology.

To characterize the predominant pathways, putative targets were searched by Kyoto Encyclopedia of Genes and Genomes analysis, which revealed the term 'chronic myeloid leukemia' (P<<0.01) with 23 target genes among the top 10 signal pathways (Table V). A disease association analysis of the predicted miRNA targets was also performed using Webgestalt demonstrating that terms as 'leukemia T-cell', 'leukemia, myeloid, acute', 'chromosome aberration', 'Precursor T-Cell Lymphoblastic Leukemia-Lymphoma', 'leukemia myeloid' 'secondary leukemia' were among the top 10 statistically significant (Table VI). In this context, the predicted target genes of miR-16-5p, miR-99b-5p and miR-144-3p were reported together with their expression values in peripheral blood leucocytes from individuals with DS and controls (Table SVIII). Expression values were obtained from a publicly available transcriptome map resulted from integration of microarray datasets through TRAM (Transcriptome Mapper) software previously (53). A total of 12 predicted target genes of miR-16-5p and 22 of miR-144-3p were identified to be differentially expressed in individuals with DS respect to controls (Ratio <0.8 and Ratio >1.2; in bold).

Table V Top 10 significant terms belonging to KEGG pathway.

Table V

Top 10 significant terms belonging to KEGG pathway.

KEGG pathway	Number of genes	%	P-value	Benjamini corrected P-value
Pathways in cancer	77	4.1	1.8×10−8	5.0×10−6
FoxO signaling pathway	35	1.9	3.6×10−7	5.0×10−5
Chronic myeloid leukemia	23	1.2	1.5×10−6	1.4×10−4
Glioma	21	1.1	4.1×10−6	2.9×10−4
Pancreatic cancer	20	1.1	1.6×10−5	9.1×10−4
Prostate cancer	24	1.3	1.7×10−5	7.8×10−4
Long-term depression	19	1.0	1.8×10−5	7.3×10−4
AMPK signaling pathway	28	1.5	8.0×10−5	2.8×10−3
cGMP-PKG signaling pathway	34	1.8	1.3×10−4	4.0×10−3
Signaling pathways regulating pluripotency of stem cells	30	1.6	1.5×10−4	4.3×10−3

[i] Bold font indicates typical Down syndrome features. KEGG, Kyoto Encclopedia of genes and genomes.

Table VI Disease association of genes targeted by miRNAs.

Table VI

Disease association of genes targeted by miRNAs.

Disease association (Term)	Number of genes	FDR
Cell transformation, neoplastic	56	7.02×10−04
Neoplastic processes	81	3.59×10−03
Leukemia, T-cell	41	3.59×10−03
Leukemia, myeloid, acute	45	3.59×10−03
Neoplasms	123	3.64×10−03
Chromosome aberrations	64	4.38×10−03
Costello syndrome	9	5.93×10−03
Precursor T-cell lymphoblastic leukemia-lymphoma	22	7.53×10−03
Leukemia, myeloid	50	9.16×10−03
Secondary leukemia	20	1.25×10−02

[i] Bold font indicates typical Down syndrome features. FDR, false discovery rate.

Furthermore, the miRNA enrichment analysis conducted by Webgestalt demonstrated that the major part (102 genes; P<<0.01) of the predicted miRNA target genes were targeted by miR-144-3p, which recognizes the binding site ATACTGT (Table VII).

Table VII miRNA enrichment analysis. Top 10 significant miRs targeting the gene list.

Table VII

miRNA enrichment analysis. Top 10 significant miRs targeting the gene list.

ID (putative binding site and miR)	Number of genes	FDR
AGCACTT, MIR-93, MIR-302A, MIR-302B, MIR-302C, MIR-302D, MIR-372,	134	0
MIR-373, MIR-520E, MIR-520A, MIR-526B, MIR-520B, MIR-520C, MIR-520D
ATACTGT, MIR-144	102	0
CAGTGTT, MIR-141, MIR-200A	124	0
GCACTTT, MIR-17-5P, MIR-20A, MIR-106A, MIR-106B, MIR-20B, MIR-519D	190	0
TGCTGCT, MIR-15A, MIR-16, MIR-15B, MIR-195, MIR-424, MIR-497	199	0
TGCACTT, MIR-519C, MIR-519B, MIR-519A	136	7.13×10−11
AAGCACT, MIR-520F	85	1.33×10−10
TTGCACT, MIR-130A, MIR-301, MIR-130B	119	8.25×10−9
ATGCTGC, MIR-103, MIR-107	75	2.31×10−8
AGCACTT, MIR-93, MIR-302A, MIR-302B, MIR-302C, MIR-302D, MIR-372, MIR-373, MIR-520E, MIR-520A, MIR-526B, MIR-520B, MIR-520C, MIR-520D	134	0

[i] Bold font indicates miRNA under study. FDR, false discovery rate; miR/miRNA, microRNA.

Among the predicted targets of miR-144-3p, two DS-associated genes, DYRK1A, encoding for dual-specificity tyrosine-(Y)-phosphorylation regulated kinase 1a and SIM1, encoding for the transcriptional factor, have been respectively reported at the 2nd and 26th positions in the Malacard gene list. DYRK1A and SIM1 gene transcripts were predicted to be targeted by miR-144-3p respectively by 8 and 6 prediction software programs (miRWalk). TargetScan (release 7.1, June 2016) predicted 2 conserved putative binding sites for DYRK1A and 1 for SIM1 gene transcripts in the 3′UTR portions.

Discussion

The present study compares the expression of miRNAs carried by biogenic nanoparticles separated from plasma samples of young individuals with DS and their healthy siblings. This is the first study that tackles noncoding RNAs encased in EVs and HDLs in association with DS pathogenesis, demonstrating specific miRNA-signatures (19). So far, only two studies have addressed an 'exosome approach' regarding this syndrome, although focusing on proteins rather than on miRNAs (54,55).

The adopted comprehensive approach of plasma NEF analysis was very informative since it allowed identification novel miRNAs probably involved in DS pathogenetic mechanisms. This approach represents a first breakthrough in respect to previous findings about miRNA and DS and sets the proof of concept for further study. It will be interesting to indicate a sharp stratification of EV and HDL miRNAs, provided the availability of sufficient starting media to perform the very challenging separation of the two nanoparticle populations.

Gene target prediction on the over-represented miRNAs identified during the initial discovery phase by microarray analysis and subsequent GO analysis of the predicted targets, suggested the possible involvement of the identified miRNAs in 'nervous system development', 'neuronal cell body' and in the pathogenesis of certain DS-associated diseases like leukemia. These are interesting results since the nervous system is remarkably affected during development and aging in individuals with DS. Furthermore, it is well known that children with DS have a markedly increased risk for a subtype of myeloid leukemia (ML) (56) classified by the World Health Organization as ML-DS (Table VIII).

Table VIII Down syndrome-associated diseases according to Malacards.

Table VIII

Down syndrome-associated diseases according to Malacards.

Related disease ID	Score
Leukemia, megakaryoblastic, with or without	12.0
Down syndrome, somatic
Acute megakaryoblastic leukemia in Down syndrome	12.0
Acute megakaryoblastic leukemia without	11.9
Down syndrome
Ayme-gripp syndrome	11.1
intellectual disability	10.4
Leukemia	10.4
Macroglossia	10.2
Dementia	10.1
Myeloid leukemia	10.1
Neuronitis	10.1

Expression profiling of the three most abundant miRNAs (miR-16-5p, miR-144-3p and miR-99b-5p) by RT-qPCR identified a trend of dysregulation in individuals with DS compared with C.

The expression of these miRNAs has not been previously associated with DS and in general little data were published regarding the expression of such miRNAs in human plasma. miR-16-5p was identified in circulating exosomes of patients with multiple myeloma (57) as well as in circulating serum exosomes of patients with esophageal adenocarcinoma (49). miR-99b-5p is known as one of five most common miRNA in human plasma exosomes by RNA deep sequencing (58).

Concerning function, miR-16-5p was previously described in rat neurons as a negative regulator of dendritic complexity and mediates brain derived neurotrophic factor (BDNF)-induced dendritogenesis by regulating the translation of the BDNF mRNA itself, supporting the hypothesis that miR-16-5p may be important during neuronal development (59).

miR-144-3p has emerged as an important miRNA implicated in certain human CNS pathologies (60,61). In a mouse model, miR-144-3p targets a number of genes implicated in the control of neuronal plasticity-associated signaling cascades, including PTEN, Spred1 and Notch1 (62).

Among the predicted targets of miR-144-3p, two DS-associated genes, DYRK1A and the transcription factor SIM1 were identified. While no experimental evidence about miR-144-3p and both DS-genes expression regulation has previously been demonstrated, it is known that DYRK1A expression can be regulated by miRNAs (63,64).

DYRK1A kinase produced by Hsa21 serves a pivotal role in the CNS and its overexpression in individuals with DS may be implicated in cognitive impairments, therefore it has been suggested as a putative therapeutic target for treating the cognitive deficiencies observed in DS (58,65,66). Deregulation of DYRK1A levels are also frequently observed in Alzheimer and Parkinson diseases (67-69).

A recent study about the transcriptome map analysis of DS vs. control human tissues (6) has demonstrated that in the differential transcriptome map of brain, induced pluripotent stem cells, blood and fibroblasts of DS and controls samples DYRK1A gene has expression ratios of 1.51, 1.36, 1.30 and 1.37, respectively. These ratios are consistent with a possible DYRK1A gene expression alteration in trisomy 21 samples due to the presence of the Hsa21 in three copies in the cells of subjects with DS. Furthermore, it is interesting to note that the SIM1 gene, which maps on human chromosome 6 is homologous to SIM2 gene mapping on Hsa21. These two genes are homologs to Drosophila single-minded (sim) gene and have expression ratios in DS vs. control brain samples of 0.75 and 0.73 respectively. It was already proposed that the human SIM1 gene may be involved in the manifestation of brain development and cognitive disability of DS (70).

Therefore, it would be interesting to experimentally verify whether DYRK1A and SIM1 genes may be targeted by miR-144-3p in order to shed light on their post-transcriptional gene regulation. Notably, the amyloid precursor protein (APP) gene, one of the most recognized genes to be overexpressed in DS (6), was not identified among the targets of the investigated miRNAs, on the basis of the selection criteria. In future studies, it would be intriguing to study if APP overexpression could be mediated by circulating miRNAs, not enclosed in nanoparticles.

The correlation of the expression levels of the three identified miRNAs with clinical features of the subjects with DS by a non-parametric test like Kruskal Wallis allowed identification of an inverse significant linear association between the development of babbling (expressed in months) and miR-16-5p and miR-144-3p in subjects with DS, which also demonstrates a high number of comorbidities. Kruskal Wallis has been applied as the present study used a small sample size where the quantitative variables analyzed, in particular miRNAs, have extremely asymmetric distributions. On the basis of the results, it can be hypothesized that miR-16-5p and miR-144-3p may be involved in controlling the timing of the babbling development: Indeed the role of miRNAs has been suggested by previous studies about speech and language impairment (71).

By means of two exploratory statistical methods, namely heatmap and cluster analysis it was demonstrated that individuals with DS expressing high levels of the three miRNAs displayed also a high number of DS-associated typical clinical features and allowed us to generate clusters of subjects well classified on the basis of the three miRNA expression data. Cluster analysis is not frequently applied to biological data, but together with heatmaps, follows the visual mining philosophy. This is an interesting strategy for understanding association and patterns into the data, visualizing simple and understandable plots.

In conclusion, the results of the present study demonstrate the importance of the comprehensive analysis of plasma miRNA vehicles, enriching for different substructures, which may lead to the discovery of unknown pathways. In the case of DS, this approach could lead to the discovery of novel molecular entities (miRNA, proteins, metabolites) that are involved in the pathology but have never been identified earlier.

Supplementary Materials

Funding

The present study was partially supported by a BIOMANE grant (University of Brescia, Brescia, Italy) to AR, AS, GDP and PB. AR was also supported from Fondo per il Finanziamento delle Attività Base di Ricerca (Ministero dell'Istruzione, dell'Università e della Ricerca) and University of Brescia research fund (ex 60%). This study was partially supported by Centro Bresciano Down (CBD), Brescia, Italy to AR and DR.

Availability of data and materials

The analyzed datasets generated during the study are available from the corresponding author on reasonable request.

Authors' contributions

AR and DR conceived the project. MC, FA, AP, CL and GC were involved in sample and medical record recruitment. AS, MV, AR wrote the manuscript. AR, LP, SB performed NEF separation and biophysical and biochemical characterization. AS, TF performed the qPCR experiments; AS, EA performed the GO analysis. GDP, GA and PB were involved in the study conception and design, analysis and interpretation of data, and critical revision of the manuscript. All authors read and approved the manuscript.

Ethics approval and consent to participate

The present study was approved by the competent Ethics Committee of Sant'Orsola-Malpighi Hospital in Bologna, approval number: 39/2013/U/Tess.

Patient consent for publication

Informed consent was obtained from the subjects for participation in this study.

Competing interests

No potential conflict of interest was reported by the authors.

Acknowledgments

Not applicable.

References