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.

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.

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.

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.

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).

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.

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.

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.

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.

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 C₁ and Cluster 2 with C₂, silhouette for the subject i contained in C₁ was: where: a(i)=average dissimilarity of I to all subjects of C₁ and b(i)=average dissimilarity of I to all subjects of C2.

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 C₁ or C₂) 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.

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.

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].

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.

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.

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).

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).

ϱ 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.

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

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 C₂ cluster, subjects with DS with high comorbidities prevailed.

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.

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).

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).

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.