Differential microRNA expression profile in blood of children with Down syndrome suggests a role in immunological dysfunction.

Down syndrome (DS), caused by trisomy of chromosome 21 (HSA21), results in a broad range of phenotypes. However, the determinants contributing to the complex and variable phenotypic expression of DS are still not fully known. Changes in microRNAs (miRNAs), short non-coding RNA molecules that regulate gene expression post-transcriptionally, have been associated with some DS phenotypes. Here, we investigated the genome-wide mature miRNA expression profile in peripheral blood mononuclear cells (PBMCs) of children with DS and controls and identified biological processes and pathways relevant to the DS pathogenesis. The expression of 754 mature miRNAs was profiled in PBMCs from six children with DS and six controls by RT-qPCR using TaqMan® Array Human MicroRNA Cards. Functions and signaling pathways analyses were performed using DIANA-miRPath v.3 and DIANA-microT-CDS software. Children with DS presented six differentially expressed miRNAs (DEmiRs): four overexpressed (miR-378a-3p, miR-130b-5p, miR-942-5p, and miR-424-3p) and two downregulated (miR-452-5p and miR-668-3p). HSA21-derived miRNAs investigated were not found to be differentially expressed between the groups. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses showed potential target genes involved in biological processes and pathways pertinent to immune response, e.g., toll-like receptors (TLRs) signaling, Hippo, and transforming growth factor β (TGF-β) signaling pathways. These results suggest that altered miRNA expression could be contributing to the well-known immunological dysfunction observed in individuals with DS.

Introduction

Down syndrome (DS) or trisomy 21 (T21) is the most common live-birth human chromosomal disorder [1]. Individuals with DS present a broad range of phenotypes including dysmorphic features [2], intellectual disability [3], congenital heart diseases [4], neurological abnormalities [5], immunodeficiency [6], and increased risk for leukemia [7], among other comorbidities.

Different hypotheses have been proposed to explain how the presence of a third copy of human chromosome 21 (HSA21) results in the DS phenotype. According to the “gene dosage effect” hypothesis, the overexpression of HSA21 genes and the downstream consequences would be directly responsible for DS characteristics [8, 9]. Besides, studies have shown the occurrence of secondary effects as a consequence of the increased gene activity on HSA21 with dysregulation of several disomic genes located throughout the genome that would affect specific cellular processes [10, 11], suggesting the existence of complex molecular mechanisms regulating RNA and protein expression. However, the determinants contributing to the complex and variable phenotypic expression of DS are still not fully known. It is likely that the phenotypic manifestations of DS result from several sources including additional content, functional elements and variability of HSA21, in addition to variability of other chromosomes throughout the genome, the chromatin structure and epigenetic modifications [12].

In light of our current understanding of the importance of genome-wide abnormalities to the phenotype of DS, we profiled the genome-wide mature microRNA (miRNA) expression in peripheral blood mononuclear cells (PBMCs) of children with DS and controls. To assess the potential biological effects of the differentially expressed miRNAs (DEmiRs) observed in our study, the potential targets of these miRNAs were predicted. In addition, following the findings of a previous study [13] by our team that compared the expression of immune-related genes between the same case and control individuals studied here, we investigated the possible interactions between those immune-related differentially expressed genes (DEGs) and the DEmiRs observed in the present study.

Materials and methods

The study protocol was approved by the Research Ethics Committee of São José do Rio Preto Medical School (CEP-FAMERP, protocol 3340/2004 and CAAE number: 88278218.4.0000.5415), in São Paulo State, Brazil. After written informed consent was obtained from parents/legal guardians, six children with DS (three males and three females, median age 4.05; range 2.30–6.10 years) and six children without DS (one male and five females, median age 3.75; range 3.60–5.40 years) were included in the study. Data on age and sex of each study participant are available in Online Resource 1.

The inclusion criteria for case and control groups were ages between 2 and 6 years, no personal history of chronic infection (such as bronchitis, asthma, or recurrent pneumonia), or actual clinical manifestations suggestive of acute infection (including the flu, cough, fever, and/or antibiotic treatment within 10 days preceding blood draw). Specific inclusion criteria were the presence of karyotypically confirmed free trisomy 21 for the case group and the absence of diseases associated with clinical manifestation of DS, such as autoimmune diseases, for the control group.

Children with DS enrolled were attending routine follow-up at the Genetics and Down syndrome Outpatient Pediatric Clinic at Hospital de Base de São José do Rio Preto (HB), in São Paulo State, Brazil. The control group consisted of children attending follow-up at the Outpatient Pediatric Clinic at HB and children attending the daycare center at the Institute of Biosciences, Letters and Exact Sciences (IBILCE), São José do Rio Preto campus—São Paulo State University (UNESP). The personal data of all participants were kept confidential.

Up to 5 mL of whole blood samples were collected by venous puncture into EDTA-containing tubes. PBMCs were isolated using Ficoll® Paque Plus (Sigma-Aldrich) density gradient centrifugation. Total RNA samples were obtained using TRIzol® Reagent (Invitrogen) following the manufacturer’s instructions. RNA quality of the twelve samples used for genome-wide mature miRNA expression analysis was assessed by the Agilent 2100 BioAnalyzer using the Eukaryote Total RNA Pico assay (Agilent Technologies) and all samples presented preserved 18S and 28S ribosomal peaks and RNA integrity number (RIN) equal to or greater than 5 [14]. RNA concentration and purity of all samples were measured by spectrophotometry using NanoDrop ND-1000 (NanoDrop Technologies).

Mature miRNA expression profiling

Mature miRNAs expression profiling was performed by Real-Time qPCR using TaqMan® Array Human MicroRNA A + B Cards Set v3.0 (Applied Biosystems), which enables accurate quantitation of 754 human miRNAs along with three reference genes (MammU6 / U6 snRNA, RNU44, and RNU48) and one negative control (ath-miR159a), according to the manufacturer’s protocol (Online Resource 2).

Briefly, total RNA (50 ng) was reverse-transcribed using the TaqMan® microRNA reverse transcription kit with Megaplex™ RT Primers Human Pools A and B (Applied Biosystems). The cDNA samples were pre-amplified using Megaplex™ PreAmp Primers, Human Pools A and B (Applied Biosystems), on a Veriti Thermal Cycler (Applied Biosystems), following the manufacturer’s instructions.

Array reactions (A and B cards) were performed in duplicate on 7900HT Fast Real-Time PCR System (Applied Biosystems) according to the manufacturer’s recommended conditions. Raw cycle quantification (Cq) data were generated by SDS RQ Manager v.1.2.1 software (Applied Biosystems) and adjusted for basal fluorescence signal and the threshold for each miRNA.

Data were filtered and Cq values above 32 were considered as undetermined, according to the TaqMan® arrays manufacturer’s instructions (Applied Biosystems). Normal distribution of data was obtained from the mean value of case and control groups separately, and Cq values outside the 75% of the central mass of the data were considered as outliers and removed from further analysis. Subsequently, when the two groups were considered together, Cq values outside the 90% of the central mass were considered outliers and excluded. Only data from genes expressed in at least four of the six biological samples in each group met our criteria to be statistically analyzed.

All expression analyses were performed in R statistical environment [15] using packages freely available under the Bioconductor project [16]. The package HTqPCR [17] was used to preprocess the data and as an interface to other packages using the qPCR data as input. miRNA expression data were normalized by geometric averaging [18] of the reference genes mmU6/U6 snRNA, RNU44, and RNU48, and relative expression (fold change, FC) was calculated by the ΔΔCq method [19] using the control group as calibrator. The identification of DEmiRs was addressed using linear models for microarray data (LIMMA) using ΔCq data as input [20]. P values were adjusted for multiple testing issues using Benjamini and Hochberg’s method [21] to control the false discovery rate (FDR) and reduce the risk of false positive. Adjusted p values less than 0.05 were considered statistically significant. Graphical data were presented as FC calculated from the mean of ΔCq of each group.

In silico functional analysis of DEmiRs

Functional characterization analysis of DEmiRs was performed using the miRNA pathway analysis web-server DIANA-miRPath v3.0 [22] (microT threshold 0.8; P value < 0.05; FDR correction; Conservative Stats), an online software suite that uses predicted miRNA targets provided by the DIANA-microT-CDS algorithm and identifies Gene Ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

The possible interactions between 17 DEGs observed in a previous study by our group [13] performed on the same individuals and the DEmiRs observed in the present study were investigated using DIANA-microT web server v5.0 [23, 24] (score threshold 0.5). Considering that miRNAs are classically known for the negative regulation of protein translation when there is a direct interaction with the target gene, only the target genes that displayed an inverse expression correlation with the miRNA were selected.

Results

Mature miRNA differential expression profiling

Of the 754 human mature miRNAs investigated using TaqMan® Arrays, 374 were detected in our experimental conditions and considered for further analysis. Initially, we observed 50 miRNAs differentially expressed between PBMCs of individuals with DS and controls (12 upregulated and 38 downregulated) (Online Resource 3). After FDR adjustment, four upregulated and two downregulated miRNAs maintained statistical significance: miR-378a-3p (FC mean = 2.53; P = 0.02), miR-130b-5p (FC mean = 2.71; P = 0.03), miR-942-5p (FC mean = 2.14; P = 0.03), miR-424-3p (FC mean = 2.36; P = 0.03), miR-452-5p (FC mean = 0.15; P = 0.04) and miR-668-3p (FC mean = 0.14; P = 0.01) (Fig. 1).

Fig. 1

Differentially expressed (DE) microRNAs (miRs) in peripheral blood mononuclear cells (PBMCs) from individuals with Down syndrome (DS) relative to non-DS children (control). The error bars represent the mean ± SEM.*FDR-adjusted P < 0.05

Eight mature forms of HSA21-derived miRNAs were investigated in the present study: miR-99a-5p and -3p, miR-155-5p and -3p, miR-125b-5p and miR-125b-2-3p, let-7c-5p and miR-802. However, only the four mature forms miR-99a-5p, miR-155-5p, miR-125b-5p, and let-7c-5p were detected in our experimental conditions. miR-155-5p was significantly overexpressed in DS in the initial analysis (Online Resource 3); however, after FDR correction, the statistical significance was lost. The other four HSA21-miRNAs did not present differential expression between DS and control groups (Fig. 2).

Fig. 2

Expression patterns of the four chromosome 21 -derived microRNAs detected in peripheral blood mononuclear cells (PBMCs) from individuals with Down syndrome (DS) relative to non-DS children (control). No statistical significance was found (FDR-adjusted P > 0.05). The error bars represent the mean ± SEM

Target genes and in silico functional analysis of DEmiRs

Table 1 presents the number of potential target genes of the six DEmiRs observed in our group of children with DS. A total of 2,401 non-redundant targets were predicted (Online Resource 4). Predicted mRNA targets are widely distributed throughout the genome and 17 of them are located on HSA21 (Table 2).

Table 1

Number of differentially expressed (DE) microRNAs (miRs) target genes predicted by the DIANA-miRPath v3.0 (DIANA-microT-CDS v5.0 algorithm)

Mature miRNA	miRBase ID Accession	Number of predicted targetsa
miR-378a-3p	MIMAT0000732	231
miR-130b-5p	MIMAT0004680	471
miR-942-5p	MIMAT0004985	1291
miR-424-3p	MIMAT0004749	110
miR-452-5p	MIMAT0001635	435
miR-668-3p	MIMAT0003881	258

aIncluding redundant targets

Table 2

Putative chromosome 21 derived-target genes of the six differentially expressed (DE) microRNAs (miRs) identified in the present study (DIANA-miRPath v3.0)

Target gene symbol	Target gene name	Target chromosome location	miRNA
BACH1	BTB domain and CNC homolog 1	21q22.11	miR-452-5p
BRWD1	Bromodomain and WD repeat domain containing 1	21q22.2	miR-452-5p; miR-130b-5p
C2CD2	C2 calcium dependent domain containing 2	21q22.3	miR-942-5p
DSCAM	DS cell adhesion molecule	21q22.2	miR-378a-3p
ERG	ETS transcription factor ERG	21q22.3	miR-942-5p
GABPA	GA binding protein transcription factor subunit alpha	21q21.3	miR-452-5p
GART	Phosphoribosylglycinamide formyltransferase, phosphoribosylglycinamide synthetase, phosphoribosylaminoimidazole synthetase	21q22.11	miR-130b-5p
HUNK	Hormonally upregulated Neu-associated kinase	21q22.1	miR-942-5p
KCNJ6	Potassium inwardly rectifying channel subfamily J member 6	21q22.1	miR-130b-5p
LCA5L	Lebercilin LCA5 like	21q22.2	miR-378a-3p
MORC3	MORC family CW-type zinc finger 3	21q22.13	miR-130b-5p
NCAM2	Neural cell adhesion molecule 2	21q21.1	miR-942-5p; miR-424-3p; miR-130b-5p
PIGP	Phosphatidylinositol glycan anchor biosynthesis class P	21q22.2	miR-942-5p
PKNOX1	PBX/knotted 1 homeobox 1	21q22.3	miR-942-5p
RRP1B	Ribosomal RNA processing 1B	21q22.3	miR-378a-3p
RUNX1	RUNX family transcription factor 1	21q22.3	miR-452-5p
SUMO3	Small ubiquitin like modifier 3	21q22.3	miR-452-5p

The in silico functional analysis including the 2401 non-redundant predicted target genes showed enrichment of genes involved in biological processes relevant to immune response (e.g., toll-like receptors—TLRs signaling pathways), coagulation, mitotic cell cycle, DNA transcription, gene expression, mRNA processing, among others (P < 0.05) (Fig. 3; Online Resource 5). KEGG analysis also showed enrichment of target genes participating in immune response pathways such a Hippo and TGF-β signaling pathways, among other pathways (P < 0.05) (Fig. 4; Online Resource 6).

Fig. 3

Gene Ontology (GO) enrichment analysis of biological processes of the differentially expressed (DE) microRNAs (miRs) targets (DIANA-miRPath v3.0; microT threshold 0.8; P < 0.05; FDR correction; Conservative Stats)

Fig. 4

Kyoto Encyclopedia of Genes and Genomes (KEGG pathway analysis of target genes of the differentially expressed (DE) microRNAs (miRs) predicted by DIANA-miRPath v3.0 (microT threshold 0.8; P < 0.05; FDR correction; Conservative Stats). The numbers within the figure refer to the number of genes involved in each KEGG pathway

The analysis of interactions between the six DEmiRs observed in the present study and the 17 immune-related DEGs previously identified [13] in the same DS and control individuals showed possible interaction between five DEmiRs and 10 DEGs (Table 3).

Table 3

Interactions between immune-related differentially expressed genes (DEGs) from our previous study [13] and the six differentially expressed (DE) microRNAs (miRs) identified in the present study (DIANA-microT web server v5.0, score threshold 0.5; TarBase v.8.0). The results are presented considering the direct effect of the miRNA on the target gene and an inverse expression correlation between them

miRNA		Target gene		Target prediction score
	Log2 FC		Log2 FC [13]
miR-452-5p	− 2.74	CCR2	1.32	0.59
miR-378a-3p	1.34	CCR7	− 1.2	0.84
		BCL2	− 0.92	0.59
miR-942-5p	1.10	IKBKB	− 0.71	0.97
		CD28	− 0.79	0.62
		CD80	− 0.94	0.62
		SKI	− 0.47	0.55
		BCL2	− 0.92	0.53
miR-130b-5p	1.44	CD28	− 0.79	0.77
		BCL2	− 0.92	0.59
		CD40	− 0.43	0.54
		IL6	− 1.74	0.53
		EDN1	− 1.79	0.53
miR-424-3p	1.24	BCL2	− 0.92	0.60

Discussion

Currently, miRBase [25, 26] version 22.1 (accessed on 08/02/2021) reports 1917 entries for stem-loop miRNA on the human genome resulting in 2654 mature miRNAs. On HSA21, there are 30 precursors reported, resulting in 38 mature miRNAs. Of these precursors, only five HSA21 sequences are classified as “high confidence miRNAs” (i.e. present high levels of certainty that they do exist) or present relevant deep sequencing data (i.e. mean of counts reads per million) (Online Resource 7). The miRNAs miR-155, miR-802, miR-125b-2, let-7c, and miR-99a, located on HSA21, have been the most extensively studied and implicated in the DS pathogenesis [27-32]. These miRNAs present two mature forms (-5p and -3p), except miR-802 [25]. In the present study, we investigated the expression pattern of all -5p and -3p mature forms of miR-155, miR-125b-5p, and miR-99a, besides miR-802; however, only -5p mature forms were detected in our samples.

The four HSA21 mature miRNAs detected in the present study did not present significant differential expression patterns between children with DS and controls after FDR adjustment. Although the overexpression of HSA21-located miRNAs has been associated with DS phenotypes [27], including low blood pressure [28], dementia [29], and altered immune processes [30, 31], there are few studies concerning miRNA expression in blood cells from individuals with DS. Xu et al. [33] identified overexpression of miR-99a, let-7c, miR-125b-2, and miR-155 in lymphocytes from children with DS by high-throughput sequencing technology and validated their data using quantitative RT-PCR. Another similar study by the same group [34] showed the same miRNAs downregulated in fetal cord blood mononuclear cells from fetuses with DS. On the other hand, comparable expression of the miRNAs miR-99a, let-7c, and miR-125b-2 were observed between cultured PBMCs from children with DS and age-matched healthy donors while overexpression of the miR-155 was observed in DS cells compared with control ones, investigated by quantitative RT-PCR [31]. Considering the different study designs, in which different methodologies for quantifying miRNAs and analyzing results were used, additional studies are necessary to better elucidate the expression pattern of HSA21 miRNAs and its consequence for the DS phenotype.

The six DEmiRs identified in the present study between children with DS and controls are not located on HSA21; four of them were observed to be upregulated (miR-378a-3p, miR-130b-5p, miR-942-5p, and miR-424-3p) and two downregulated (miR-452-5p and miR-668-3p). Although some DEmiRs identified in the present study have been previously found to be dysregulated in other DS tissues [35-37], to our knowledge, this is the first study to report altered expression of these six DEmiRs in PBMCs from individuals with DS. Altered expression of miRNAs not located on HSA21 with possible involvement in DS phenotypes has been previously reported [33-40], supporting the existing hypothesis of secondary transcriptional changes as a consequence of the trisomy 21 [41].

According to the literature, sex can influence gene expression levels [42, 43], as well as the expression of miRNAs [44]. Considering that in the present study the sex was not matched between individuals with DS and controls, we investigated if there was difference in the expression of the six DEmiRs and the four HSA21 mature miRNAs between male and female participants regardless the study group (case or control). miR-452-5p and miR-668 were observed to be overexpressed in females (data not shown), indicating a possible sex bias in miRNA expression. Therefore, we compared the expression of these two miRNAs between case and control groups considering only female participants (data not shown) and obtained the same results from the total sample (both downregulated in the DS group), indicating that there is no sex bias in our results. Data on sex and age for each sample include in the study are available in the Online Resource 1.

Enrichment analyses of GO biological processes and KEGG pathways including the 2401 non-redundant predicted targets of the six DEmiRs point to the contribution of these genes to pathways related to the immune system, including the TLRs, TGF-β, and Hippo signaling pathways. TLRs are key elements of innate immunity as they participate in the recognition of pathogens, antigen presentation, apoptosis, and production of interferons by the virus-infected cells [45]. TLR signaling is dysregulated in children with T21 and this impairment is believed to contribute to the chronic inflammation and sepsis observed in DS [46]. The TGF-β signaling pathway is a pivotal regulator of immune responses, participating in the proliferation, survival, activation, differentiation, and repertoire diversity of B- and T- cells under normal conditions and immune challenges [47]. Moreover, TGF-β regulates the development and functionality of innate cells (dendritic cells, macrophages, granulocytes, and natural killer cells) and the peripheral tolerance against self-antigens [47]. In DS, increased TGF-β levels have been previously associated with transient abnormal myelopoiesis and complications of this condition [48] and pointed as a candidate biomarker for prenatal diagnosis [49].

Studies have revealed extensive roles for the Hippo signaling pathway in adaptive and innate immune systems regulation. Components of the Hippo pathway (MST1/2) are involved in T-cell development and differentiation, B cell homeostasis, and they play key roles in antigen-presenting cells, such as macrophages and dendritic cells [50]. Also, there is evidence of crosstalk between Hippo and TLR signaling pathways [51], along with other signaling networks involved in immune regulation [52]. To our knowledge, there are no studies concerning the role of the Hippo pathway in the DS phenotype and our study is the first to point out a possible association between Hippo pathway dysregulation and DS.

Several of the predicted HSA21-located target genes of our DEmiRs play a role in the immune system. BACH1 is involved in the regulation of the expression of macrophage-associated genes, B-cell development, and antigen presentation [53]; BRWD1 cooperates with networks to drive B-cell development [54]; ERG participates in the control of B lymphopoiesis [55]; GABPA is a regulator of cytokine secretion [56]; RUNX1 regulates TLR signaling pathways and inflammatory cytokine production [57]. Besides the fact that they are present in triplicate in individuals with DS, the altered expression of miRNAs that target these genes could be modulating their expression pattern giving further support to the role of these miRNAs in DS pathogenesis.

Previous studies from our research group showed dysregulation of immune and inflammation-related genes non-located on HSA21 in children with DS [13, 58]. The study of immune-related genes [13] was carried out in the same cohort of individuals with DS and controls included in the present study; thus, we deemed it interesting to investigate a possible interaction between the immune-related DEGs previously identified and the six DEmiRs identified in the present study. Fourteen possible miRNA–target interactions were found (Table 3). These immune-related genes are involved in key immune response processes, such as regulation of T and B cells, production of cytokines, antibodies and autoantibodies, activation of signaling pathways, and their dysregulation could negatively influence the immune response in DS, as previously discussed [13].

Immune dysregulation has been widely described in DS [6, 59–61]. In line with the contribution of immunological alteration and infectious diseases to the pathophysiology of DS, individuals with T21 present increased morbidity and mortality due to respiratory tract infections [62, 63]. Since the start of the COVID-19 pandemic caused by SARS-CoV-2, there is concern about the individuals with DS being a risk group for the disease [64]. In fact, studies have shown that they present an increased risk for severe disease course and COVID-19–related hospitalization and death [65-67]. Recent network analysis of DS and SARS-CoV-2 revealed that the susceptibility to infection, worse prognosis, and complications of COVID-19 could be related to genetic factors associated with DS that predispose to unfavorable conditions, such as the dysregulation of genes involved in the immune response [68].

Our findings suggest that altered miRNAs expression patterns could be contributing to the immunological dysfunction observed in children with DS by targeting genes involved in several key processes crucial for the immune system to function properly. We emphasize the importance of investigating and understanding the immune system of DS individuals who constitute a clinically vulnerable population so that they can have better support and clinical management, and consequently  be better  in a future pandemic scenario. Considering that one coding gene can be regulated by multiple miRNAs [69] and one miRNA can regulate multiple targets [70], studies performing high throughput screening of miRNA and mRNA expression patterns could contribute to the understanding of the gene expression regulation in DS. Possible miRNA–target interactions identified in the present study are candidates for further investigation and experimental validation to understand the impact of DEmiRs on the mRNA expression of specific genes and consequently on the pathophysiology of DS.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 13 KB)

Supplementary file2 (XLSX 85 KB)

Supplementary file3 (PDF 214 KB)

Supplementary file4 (XLSX 143 KB)

Supplementary file5 (XLSX 14 KB)

Supplementary file6 (XLSX 12 KB)

Supplementary file7 (XLSX 15 KB)