This article is mentioned in:

Introduction

Down syndrome is one of the most frequently occurring chromosomal abnormalities in humans, with a worldwide incidence of ~1:1,000 births annually (1). This disease is caused by the presence of part or a full extra copy of chromosome 21 and therefore is also called trisomy 21 (2). Down syndrome is typically associated with craniofacial abnormalities, intellectual disabilities, physical growth delays and a variety of distinctive physical features, including slanted eyes and a large tongue (3). Furthermore, sufferers of Down syndrome have an increased susceptibility to certain diseases, including leukemia, Alzheimer's disease, congenital heart disease and thyroid disorders (4). Nowadays, bioinformatics analysis has been widely performed to investigate disease-associated genes and genes on human chromosome 21, as well as other chromosomes that have been linked to Down syndrome (5,6). However, the roles of microRNAs (miRs) and miR-mRNA connections in Down syndrome are largely unknown.

miRs represent a class of single-stranded noncoding RNAs, which serve an important regulatory role in a variety of biological processes including differentiation, proliferation and apoptosis, through post-transcriptionally regulating gene expression via binding with complementary targets (7,8). By performing bioinformatic analyses, several known and novel miRs that are differentially expressed in Down syndrome have been identified. For example, miR-1246 has been identified as a link between the p53 family and Down syndrome (9), and miR-138-5p and its target zeste homolog 2 have been demonstrated to function in the hippocampus of Down syndrome patients (10).

The extra genetic material on chromosome 21 in Down syndrome subjects, in particular the material involving a critical region on chromosome 21, leads to a global genomic dysregulation of all chromosomes, including abnormalities of thymic structure and function (11). The immune system in Down syndrome is intrinsically deficient. Down syndrome patients present abnormal thymuses characterized by lymphocyte depletion, cortical atrophy and loss of corticomedullary delimitation (12,13). To investigate epigenetic mechanisms associated with the thymic gene co-expression network and Down syndrome genomic dysregulation, Moreira-Filho et al (14) performed miR target analyses for differentially expressed miRs using a network-based approach. RNA sequence (RNA-seq) data, including miR (GSE69210) and mRNA (GSE70573) data collected from the thymic tissues of both Down syndrome and karyotypically normal subjects, were deposited in a public repository. In the present study, miR and mRNA databases were downloaded to establish a comprehensive association between miR expression and target mRNA expression in thymus tissues in Down syndrome, as an initial step towards further understanding the regulation of gene expression in this disease.

Materials and methods

Data collection and processing

The mRNA expression profile of GSE69210 (14), containing 10 thymus samples collected from Down syndrome subjects and 10 thymus samples from karyotypically normal subjects, was downloaded from the public functional genomics data repository Gene Expression Omnibus (GEO) of the National Center of Biotechnology Information (www.ncbi.nlm.nih.gov/geo/), platform GPL13497. The miR expression profile of GSE70573 (14) was downloaded from the GEO database, platform GPL13263 and contained eight human thymus samples from Down syndrome subjects, and eight thymus samples from karyotypically normal subjects.

According to the names of the thymus samples and the clinical information of subjects, paired miR and mRNA expression profiles were obtained from eight Down syndrome subjects and eight karyotypically normal subjects. The raw probe-level data of both miR and mRNA were converted into expression measures. The expression values of all probes for a given gene were reduced to a single value by taking their average expression value. Then, RNA-seq data were log2-transformed and median normalized using the limma package in R (version 3.3.3) (15,16).

Weighted gene co-expression network analysis (WGCNA)

Following data normalization, the miR and mRNA modules associated with Down syndrome using WGCNA in R were screened (15). First, a measure of similarity between miR or mRNA expression profiles that measures the level of concordance between expression profiles was defined. The absolute value of the Pearson correlation for each pair of genes was used to denote this similarity measure. Then, the similarity matrix was transformed into an adjacency matrix following definition of an adjacency function based on scale-free topology criterion. The Pearson correlation was transformed into a matrix of the power adjacency functions (so-called soft threshold) using the following formula: Power=|correlation|β. Next, the parameters indicating the adjacency function, including the sign function and the power function, were determined. Finally, modules associated with Down syndrome were identified according to the topological overlap matrix-based dissimilarity measure in conjunction with a clustering method.

Identification of differentially expressed genes (DEGs) and miRs (DEmiRs)

Differential miR or mRNA expression between Down syndrome samples and karyotypically normal samples was determined from the normalized data using the limma package for R, with P<0.05 and |log2 Fold Chance (FC)|>0.585 set as thresholds. The identified DEmiRs and DEGs were mapped to the miR and mRNA modules screened by WGCNA, and overlapped miRs and mRNAs were defined as DEGs or DEmiRs.

Bidirectional hierarchical clustering analysis

Bidirectional hierarchical clustering (16) of DEGs or DEmiRs was performed based on normalized expression values by employing the Euclidean distance method (17) and the heatmap package in R (18). The results were then visualized on a heatmap. Bidirectional clustering includes expression profile clustering as well as sample clustering. Expression profile clustering can be used to link genes with similar expression profiles together, which can aid examination of the data; sample clustering based on expression levels can help visually identify whether screened DEGs and DEmiRs have sample specificity.

Co-expression network construction of DEmiR-targeted gene pairs

DEmiR target sites were retrieved using a public miRanda algorithm (19) found at www.microRNA.org. Following target site retrieval, target mRNAs were mapped to the mRNA modules screened by WGCNA to identify DEmiR-targeted gene pairs.

Then, the Pearson correlation coefficients of DEmiR-targeted gene pairs were calculated based on the expression levels using the Cor function in R. Gene pairs with correlation coefficient >0.6 were kept. Following this, the co-expression network of DEmiRs and screened DEmiR-targeted gene pairs was constructed and visualized using Cytoscape (version 3.4.0) (20). The modules in the co-expression network were identified and annotated using MCODE (degree cutoff=2, node score cutoff=0.2 and k-core=2) and BiNGO plugins (adjusted P-value<0.05) (21,22).

Gene Oncology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment

The functional characteristics of the genes in the co-expression network were identified by GO analysis (gostat.wehi.edu.au/) (23). During GO analysis, genes were connected to biological GO annotations. The most enriched GO terms were statistically identified from thousands of GO annotation terms by enrichment analysis based on the hypergeometric distribution algorithm as follows, with P<0.05.

N represents the number of genes with GO-based functional annotation; K represents the number of DEGs out of genes with GO-based functional annotation; M represents the number of genes that were annotated to the specific GO term; i represents the specific GO term.

Additionally, KEGG orthology-based annotation system (kobas.cbi.pku.edu.cn/home.do) (24) analysis was applied to predict the enriched pathway in the KEGG database of genes in the co-expression network, based on the cumulative hypergeometric distribution algorithm with P<0.05.

Results

Identification of Down syndrome-associated miRs and mRNAs based on WGCNA

Following data normalization, the paired miR and mRNA expression profiles collected from eight Down syndrome subjects and eight karyotypically normal subjects were analyzed using WGCNA.

Definition of network adjacency function

By applying approximate scale-free topology, the value of power (also termed the soft threshold), a weighted parameter of the network adjacency matrix, was investigated. The definition of threshold parameters for network construction and calculation of scale-free degree distributions is demonstrated (Fig. 1). Higher values of the square of correlation coefficient indicated a closer approximation to scale-free distribution of the network. The value of power was chosen to be the value of the square that first reached 0.9, which was 10 and 4 for miR and mRNA data, respectively, as illustrated in Fig. 1A and B.

Figure 1. Screening for (A) miRNAs and (B) mRNAs associated with Down syndrome using weighted gene co-expression network analysis. An investigation of the value of power, a weighted parameter of the network adjacency for (Aa) miRNAs and (Ba) mRNAs. The horizontal axis represents the power and the vertical axis represents the square of the correlation coefficient; the red line indicates the square of the correlation coefficient of 0.9. Dendrograms of screened (Ab) miRNA or (Bb) mRNA modules associated with Down syndrome. The different colors in the bottom row represent the different modules. The correlation of the modules with Down syndrome for (Ac) miRNAs and (Bc) mRNAs. The horizontal axis represents different modules in different colors; the vertical axis represents the correlation coefficient between genes in modules and Down syndrome; the horizontal line represents the correlation coefficient of 0.8. miRNAs, microRNAs.

Network construction and module definition

After the network was established, the topological overlap matrix was converted to a dissimilarity measure and submitted to hierarchical clustering, creating a dendrogram that clusters similarly expressed genes into discrete branches. Then, modules containing at least 30 miRs and 100 mRNAs were assigned using mixed dynamic tree-cutting algorithm criteria. Next, height cutoffs of 0.95 were selected for miR and 0.90 for mRNA to cut branches off the tree. A total of 6 miR modules and 20 mRNA modules were obtained. The resulting branches corresponding to the miR and mRNA modules are presented in Fig. 1Ab and Bb, and replaced by the bottom row color-code.

Identification of Down syndrome-associated miR and mRNA modules

The correlation coefficient between the eigenvalue of each color module and the gene difference between Down syndrome samples and normal samples, as well as the significance of correlation, were calculated. As demonstrated in Fig. 1Ac and Bc, the modules with correlation coefficients >0.8 were selected. For miR modules, only one module (blue) containing 76 miRs was obtained; for mRNA modules, a total of eight modules (black, cyan, grey, light yellow, magenta, purple, red and salmon) containing 888 mRNAs were selected.

Identification of DEmiRs and DEGs

For paired miR and mRNA expression profiles, a total of 428 DEGs and 54 DEmiRs were obtained, at P<0.05 and |log2 FC|>0.585. Following mapping of these miRs and mRNAs to the miR and mRNA modules screened using WGNCA, 12 downregulated DEmiRs and 237 DEGs (26 downregulated and 211 upregulated) were obtained, all of which were significantly associated with Down syndrome (P<0.05; data not shown).

Bidirectional hierarchical clustering analysis of DEGs and DEmiRs

According to gene expression patterns, 237 DEGs and 12 DEmiRs (downregulated) were clustered through bidirectional hierarchical clustering analysis. Then, the Down syndrome samples could be distinguished clearly from the control samples (Fig. 2A for miR and Fig. 2B for mRNA).

Figure 2. Bidirectional hierarchical clustering and heat map of differentially expressed (A) miRs and (B) genes. miR, microRNA; DS, Down syndrome; CT, control samples.

DEmiR-targeted gene co-expression network

The target sites of 12 DEmiRs were retrieved from the miRanda database and blasted with 237 DEGs. Then, 255 DEmiRs-DEGs pairs, including 6 DEmiRs and 106 DEGs were obtained. Next, the expression correlation between two genes out of the 106 DEGs was evaluated at the expression correlation coefficient >0.9; a total of 231 gene pairs were obtained. Finally, a DEmiR-regulated gene co-expression network was constructed. As demonstrated in Fig. 3, there were 112 nodes containing 6 downregulated miRs, 9 downregulated mRNAs and 97 upregulated mRNAs, and 486 edges (255 miR-gene connections and 231 gene-gene connections).

Figure 3. Gene co-expression network regulated by DEmiRs. Square nodes represent miR; the regular triangles and downward facing arrowheads represent the upregulated and downregulated mRNAs, respectively; green lines represent connections between miRs and mRNAs; black and red lines represent connections between genes positively and negatively correlated, respectively; grey nodes represent genes involved in functional modules. DEmiRs, differentially expressed microRNAs.

Functional enrichment analysis

GO and KEGG pathway enrichment analyses were performed for genes in the DEmiR-regulated gene co-expression network. Genes were enriched in multiple GO categories including microsome, vesicular fraction and actin binding, and 11 KEGG pathways, as demonstrated in Fig. 4.

Figure 4. Functional enrichment analysis. (A) Functional annotations of Gene Ontology terms. (B) Kyoto Encyclopedia of Genes and Genomes pathways enrichment. The values in brackets represent P-values.

Functional division of DEmiR-regulated gene co-expression network

The co-expression network was further divided into several functional modules. A total of two DEmiR (miR-145 and miR-183)-regulated modules were then obtained (Figs. 5 and 6). The genes in these two modules, which were enriched in nine GO categories including cellular lipid metabolic process, lipid metabolic process and cellular response to stimulus, and five GO categories including response to stimulus or stress were functionally annotated, as illustrated in Tables I and II.

Figure 5. miR-145-regulated functional modules. Square nodes represent miRs; regular triangles represent upregulated mRNAs; dotted lines represent connections between miRs and mRNAs; solid lines represent connections between positively correlated genes. miRs, microRNAs. FCER1A, high affinity immunoglobulin ε receptor subunit α; DOK7, downstream of tyrosine kinase 7; FLG2, filaggrin-2; PSAPL1, proactivator polypeptide-like 1; DLEC1, deleted in lung and esophageal cancer protein 1; SHROOM3, shroom family member 3.

Figure 6. miR-183-regulated functional modules. Square nodes represent miR; regular triangles and downward facing arrowheads represent upregulated and downregulated mRNAs, respectively; dotted lines represent connections between miRs and mRNAs; solid lines represent connections between positively correlated genes. miR, microRNA. CHAF1B, chromatin assembly factor 1 subunit B; HSPA13, heat shock protein family A member 13; ROBO2, roundabout guidance receptor 2; GPX8, glutathione peroxidase 8; P2RY2, purinergic receptor P2Y2; REEP3, receptor accessory protein 3; XRRA1, X-ray radiation resistance associated 1.

Table I. Gene Ontology annotation of functional modules regulated by microRNA-145.

Table I.

Gene Ontology annotation of functional modules regulated by microRNA-145.

GO-ID	Description	P-value	Genes in test set
GO:0044255	Cellular lipid metabolic process	8.41×10−3	PSAPL1, FCER1A
GO:0006629	Lipid metabolic process	1.88×10−2	PSAPL1, FCER1A
GO:0051716	Cellular response to stimulus	2.60×10−2	SHROOM3, FCER1A
GO:0051179	Localization	1.94×10−2	SHROOM3, FCER1A
GO:0050794	Regulation of cellular process	2.22×10−2	SHROOM3, FCER1A, DLEC1
GO:0050789	Regulation of biological process	2.52×10−2	SHROOM3, FCER1A, DLEC1
GO:0050896	Response to stimulus	2.68×10−2	SHROOM3, FCER1A
GO:0065007	Biological regulation	2.91×10−2	SHROOM3, FCER1A, DLEC1
4GO:004237	Cellular metabolic process	4.35×10−2	PSAPL1, FCER1A

[i] FCER1A, high affinity immunoglobulin ε receptor subunit α; DOK7, downstream of tyrosine kinase 7; FLG2, filaggrin-2; PSAPL1, proactivator polypeptide-like 1; DLEC1, deleted in lung and esophageal cancer protein 1; SHROOM3, shroom family member 3.

Table II. Gene Ontology annotation of functional modules regulated by microRNA-183.

Table II.

Gene Ontology annotation of functional modules regulated by microRNA-183.

GO-ID	Description	P-value	Genes in test set
GO:0050896	Response to stimulus	4.99×10−3	ROBO2, CHAF1B, P2RY2, XRRA1, GPX8
GO:0051716	Cellular response to stimulus	5.59×10−3	ROBO2, CHAF1B, P2RY2
GO:0042221	Response to chemical stimulus	1.69×10−2	ROBO2, P2RY2, GPX8
GO:0006950	Response to stress	2.84×10−2	CHAF1B, P2RY2, GPX8
GO:0023052	Signaling	1.23×10−2	ROBO2, P2RY2, D4S234E

[i] CHAF1B, chromatin assembly factor 1 subunit B; HSPA13, heat shock protein family A member 13; ROBO2, roundabout guidance receptor 2; GPX8, glutathione peroxidase 8; P2RY2, purinergic receptor P2Y2; REEP3, receptor accessory protein 3; XRRA1, X-ray radiation resistance associated 1.

Discussion

Down syndrome, or trisomy 21, is one of the most frequently occurring chromosomal abnormalities in humans (2). To better understand the regulatory mechanism of gene expression in Down syndrome, miR target gene expression profiles were systematically analyzed in the present study.

miR and mRNA expression profiles in the thymus tissue of human Down syndrome patients were first downloaded from GEO, a public functional genomics data repository, from a study by Moreira-Filho et al (14). These RNA-seq data were preliminarily screened using WGCNA. WGCNA is a powerful statistical method of constructing co-expression networks describing correlation patterns among genes across microarray samples (14). In the present study, WGCNA, as implemented in R, was used to identify Down syndrome-associated miR and mRNA modules. Notably, DEGs and DEmiRs based on the expression levels of miRs and mRNAs between Down syndrome thymus samples and normal control samples were also identified. Then, these DEGs and DEmiRs were mapped to the screened miR or mRNA modules by WGCNA and were identified as DEGs and DEmiRs. Next, a gene co-expression network regulated by six DEmiRs, including miR-30c, miR-139-5p, miR-183, miR-145, miR-150 and miR-135b, was constructed. Additionally, DEmiR target genes were enriched in various biological processes and signaling pathways, including those modulating response to stimulus and response to stress.

miRs are known to post-transcriptionally regulate gene expression by binding to complementary targets which are involved in a wide range of important biological processes. They also can be used as predictive tools and intervention targets for several diseases. Previous studies have been performed investigating the miRs and their predicted mRNAs in Down syndrome. Using bioinformatics tools, Lim et al (25) profiled genome-wide miR expression in placenta samples from euploid or Down syndrome fetuses and predicted the functions of DEmiRs. A total of 34 placental miRs were identified that were differentially expressed in euploid fetuses with Down syndrome and proposed that a number of Down syndrome-related genes on chromosome 21 are predicted targets of these miRs. In the study of Li et al (26), the miRs from the brain and peripheral tissues of Downs patients were screened and demonstrated an upregulation of miR-155 and a decrease in the abundance of the miR-155 mRNA target complement factor H (CFH). Therefore, it was hypothesized that a generalized miR-155-mediated downregulation of CFH may contribute to the immunopathological deficits associated with Down's syndrome. In the hippocampus of the Ts65Dn Down syndrome mouse model, the genome-wide regulatory effects of miR-155 and miR-802, which are overexpressed in Down syndrome individuals, have been demonstrated by the lentiviral miR-sponge strategy (27). An integrated analysis of miR and mRNA expression in the hippocampus has been presented by Shi et al (10). They identified a number of DEmiRs and their predicted mRNAs. Of these identified miR-mRNA interactions, miR-138 and its target histone-lysine N-methyltransferase EZH2 have been demonstrated to function in the hippocampus in Down syndrome (10). In addition, Moreira-Filho et al (14) used whole thymic tissue from Down syndrome, karyotypically normal subjects, and a network-based approach for gene co-expression network to identify modular transcriptional repertoires and the interactions between all of the system's constituents. Down syndrome is described in the study as a thymic tissue adaptation under trisomy 21 genomic dysregulation, and this adaption may be driven by epigenetic mechanisms acting at the chromatin level and through miR-based control of transcriptional programs (14).

In the present study, the DEmiR and DEmiR-regulated gene co-expression network in the thymus tissue of Down syndrome patients suggested an association between these identified miRs/mRNAs and the pathogenesis of Down syndrome, as well as the development of several complications in Down syndrome patients. A previous study identified that through binding to the 3′untranslated region of miR-30c, miR-30c can bind and post-transcriptionally attenuate the expression of runt-related transcription factor 1 (RUNX1), a key regulator of definitive hematopoiesis in the embryo and adults, which resides on chromosome 21 (28). Due to the increased gene dosage of RUNX1, patients with Down syndrome have an increased risk of developing acute megakaryoblastic leukemia (29). In addition, patients with Down syndrome have decreased expression of miR-30c, and increased expression of RUNX1. In the present study, miR-30c was demonstrated to target several other sites linked to Down syndrome in addition to RUNX1, including heat shock 70 kDa protein 1B (HSPA1B), protein BTG3 (BTG3), double-stranded RNA-specific editase 1 (ADARB1), pericentrin (PCNT), DENN domain-containing protein 5B (DENND5B), α-2-macroglobulin-like protein 1 (A2ML1) and ectoderm-neural cortex protein 1 (ENC1). Patients who have Down syndrome are at high risk of developing early-onset Alzheimer's disease. HSPA1B, has been demonstrated to be associated with Alzheimer's disease pathophysiology (30). Additionally, as reported previously, there is an increased expression of heat-shock proteins in patients with both Alzheimer's disease and Down syndrome (31). These results suggest that HSPA1B may be involved in the onset of Alzheimer's disease in Down syndrome patients. BTG3, a candidate tumor suppressor gene induced by p53 following DNA damage, which is also located on chromosome 21 is significantly altered in adult Down syndrome patients (32). As previously described, BTG3 appears to be tightly regulated in Down syndrome (2). ADARB1 is known to serve a key role in brain development and function, and has been reported to demonstrate a dosage imbalance in Down syndrome mouse models at the transcriptional level (33). Furthermore, previous RNA expression analysis identified that PCNT on chromosome 21, which is involved in Down syndrome, is dose-sensitive (34). Furthermore, DENND5B, A2ML1 and ENC1 have been identified as differentially expressed in fibroblasts carrying trisomy 21 compared with normal fibroblasts, following transcriptome analysis of induced pluripotent stem cells from monozygotic twins (35). The results of the present study, together with evidence from previous research, indicates that miR-30c-regulated gene co-expression networks may serve an important role in the development of detrimental phenotypes in Down syndrome patients, particularly regarding Alzheimer's disease.

In the present study, the gene co-expression network was further divided and miR-145, and miR-183 regulated functional modules, and genes involved in these two modules, were demonstrated to be associated with responses to stimuli and responses to stress. miR-145 is known as a tumor-suppressive miR that is downregulated in various types of cancer, including prostate, bladder, colon and ovarian cancers (36). miR-183 is located on human chromosome 7 and has been implicated in several key cellular functions including neurosensory development (37). It is also reported to function as an oncogene by targeting early growth response protein 1 and promoting tumor cell migration (37). A previous study demonstrated that oxidative stress occurs in the pathogenesis of Alzheimer's disease with Down syndrome due to a deregulation of gene/protein expression associated with trisomy 21 (38). Increased production of reactive oxygen species accompanied by mitochondrial dysfunction can trigger Alzheimer's disease onset in Down syndrome patients (39). In addition, it has been reported that levels of oxidative stress and induction of heat-shock protein response in the amniotic fluid were higher in women carrying Down syndrome fetuses compared with non-Down syndrome controls (38). Additionally, P2Y purinoceptor 2, regulated by miR-183, and widely expressed by thymic epithelial cells (40), was identified by Moreira-Filho et al (14) following miR target interaction analysis in the thymus of Down syndrome patients.

Although the present study presented an integrated analysis of the miR-mRNA network in Down syndrome using bioinformatics tools, the study is faced with an inherent limitation; the lack of biological validation of the bioinformatics results. The results obtained from this bioinformatics analysis should be confirmed in further experimental studies.

In conclusion, the present study represents an integrated analysis of miR and mRNA expression in the thymus in Down syndrome patients. A gene co-expression network of six miRs and 122 target mRNAs was constructed. These miRs and their targets, particularly those of miR-30c, miR-145 an miR-183, may be involved in the pathogenesis and development of complications in Down syndrome. However, these results should be confirmed in further experimental studies.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

MC and LS participated in the design of the present study and performed the statistical analysis. XHao, MZ, LZ and JB conducted the study and analyses, and collected important background information. XHan and CG drafted the manuscript, and analyzed and interpreted the data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References