Temporal dynamics of miRNA expression during cerebral cortex developmentTo study the expression patterns of miRNAs during mammalian corticogenesis, we performed small bulk RNA sequencing in both female and male mice at a progenitor-dominated developmental phase (embryonic day E14), at an intermediate time point (E17) and at birth (postnatal day P0) when neurons from all six layers of the neocortex have been born27,28 (Fig. 1). Using stringent filtering criteria for expression levels (>10 CPM in at least five samples), we detected a total of 489 miRNAs with only 13 of those corresponding to novel miRNAs. In agreement with previous reports17,29, miR-9, members of the let-7 family, miR-128, and miR-124 were among the top 20 most abundantly expressed miRNAs in the developing cerebral cortex (Supplementary Fig. 1). Principal component analysis (PCA) revealed a distinct separation of samples according to their developmental stage. However, E17 and P0 samples had a more similar global expression pattern to each other compared to E14 (Fig. 2a). Since we did not observe apparent differences related to sex, female and male samples were processed together in the subsequent analyses.Fig. 1: Schematic representation of the study.To investigate the expression patterns of miRNAs during cerebral cortical development in the mouse, we performed RNA-seq of bulk tissue at E14, E17 and P0 as well as in NPCs isolated from the cerebral cortex and differentiated into neurons. Subsequently, we employed weighted gene co-expression network analysis (WGCNA) to identify modules of co-expressed miRNAs. We then constructed a network of miRNAs sharing more target genes than expected by chance. Co-targeting relationships between selected miRNAs were validated using luciferase reporter assays, demonstrating enhanced gene silencing effect on neurodevelopmentally-relevant genes for combinations of miRNAs compared to individual miRNAs.Fig. 2: Expression changes of miRNAs during mouse cerebral cortex development.a Principal component analysis of bulk small RNA sequencing samples from the cerebral cortex of female and male mice at different developmental stages. n = 4 biological replicates for females at E14, n = 6 in all remaining groups. b Venn-diagram of upregulated miRNAs at earlier compared to later stages of cortical development. c Venn-diagram of downregulated miRNAs at earlier compared to later developmental stages. d–f Volcano plots of differentially expressed miRNAs. Positive log2 fold changes indicate miRNAs upregulated at earlier compared to later developmental stages in each pairwise comparison (E14 vs. E17, E14 vs. P0 and E17 vs. P0). g Heatmap of miRNAs differentially expressed in each pairwise comparison of different developmental stages. Mean expression values per developmental stage are shown as z-scores.We detected the highest number of differentially expressed miRNAs at the E14 vs. P0 time points, where 169 miRNAs were up- and 172 were down-regulated (Fig. 2b–f, Supplementary data 1). However, the major transcriptional shift already occurred in the transition from E14 to E17, where ~57% of all detected miRNAs significantly changed their expression level (122 up- and 144 down-regulated miRNAs). In contrast, only 36% of all detected miRNAs were differentially expressed in E17 vs. P0 samples (82 upregulated and 93 downregulated miRNAs). Interestingly, a higher number of miRNAs were consistently downregulated at earlier developmental time-points in all pairwise comparisons (45 vs. 32 upregulated miRNAs, Fig. 2b, c, g).To study the expression of miRNAs specifically during neurogenesis, we isolated NPCs from the cortices of mouse embryos at E14 and differentiated them in vitro into neurons. Small RNA sequencing of these cells revealed 110 miRNAs that were upregulated and 113 miRNAs that were downregulated in NPCs compared to neurons, respectively. miR-124-3p and miR-369-3p were the most significantly upregulated miRNA in neurons, whereas miR-155-3p and miR-34b-3p were the most significantly upregulated miRNAs in NPCs (Supplementary Fig. 2a, d, Supplementary data 1). 66% of the upregulated (72 out of 110) and 65% of the downregulated (73 out of 113) miRNAs in NPCs vs. neurons were also differentially expressed in the comparison of bulk samples from E14 vs. P0, confirming the replicability of the expression patterns we observed (Supplementary Fig. 2b–d). Interestingly, these overlapping miRNAs corresponded to 43% of all upregulated and 42% of the downregulated miRNAs in E14 vs. P0, suggesting that a large proportion of the differentially regulated miRNAs might be attributable to the non-neuronal cells in the bulk samples. This observation is plausible, as gliogenesis takes place roughly between E17 and P0 of cortical development28,30 therefore the bulk samples from P0 likely contained both neuronal and non-neuronal cell types. In agreement, several of the developmental-stage-dependent expression changes we observed align with previously reported roles of miRNAs in inhibiting or promoting gliogenesis. For instance, miR-106a-5p was upregulated at E14 compared to E17 and P0 (Fig. 2g). miR-106a was previously reported to suppress gliogenic differentiation of neural stem cells and induce neurogenic cell fate commitment31, which fits well with the upregulation of expression during peak neurogenesis in our study. Furthermore, we observed an upregulation of the gliogenic miR-33832 specifically at P0 but not in neurons compared to NPCs (Supplementary data 1), thus underscoring the potential role of this miRNA specifically in promoting differentiation of glia cells. As another example, miR-153 whose inhibition was shown to confer gliogenic competence to neural stem cells33, was upregulated in neurons compared to NPCs in our analysis (Supplementary data 1).Weighted gene co-expression gene network analysis identifies sets of co-regulated miRNAs during the development of the cerebral cortexNext, we applied weighted gene co-expression network analysis (WGCNA) to construct networks of co-expressed miRNAs and better characterize their transcriptional dynamics along the cortical developmental trajectory (Fig. 1). We obtained 18 modules that were assigned to arbitrary colors (Fig. 3a, detailed information of the specific miRNAs included in each WGCNA module is included in Supplementary data 2). Module eigengenes (MEs) that correspond to the first dimension in a principal component analysis of the expression matrix of the corresponding WGCNA module, served as proxies for the characteristic transcriptomic signature of each module. Six MEs were significantly negatively correlated with the developmental time point indicating that the miRNAs in the respective modules showed overall reduced expression levels at later compared to earlier developmental stages (Fig. 3b). In contrast, four modules contained miRNAs whose eigengene expression significantly increased during brain development. The black module, which contained the highest number of miRNAs (139), had the strongest negative correlation with the developmental stage. The overall expression pattern captured by the ME, indicated that the black module contained miRNAs whose expression level peaked at E14 and dropped at E17, remaining stable afterwards (Fig. 3c). The green module was the second largest, containing 101 miRNAs with generally linearly increasing expression levels from E14 to P0 (Fig. 3d). Remarkably, although these two modules were regulated in opposite directions, they shared most of their targets (6244 common targets, corresponding to 79% of all black module and 72% of all green module targets, Supplementary Fig. 3c). Gene ontology analysis using the genes targeted by the miRNAs revealed that both the black and the green modules are involved in the regulation of key nervous system developmental processes such as axon development and guidance, neuron projection, dendrite development, and synapse organization (Fig. 3e, f, Supplementary data 3).Fig. 3: Weighted miRNA co-expression analysis.The co-expression network was derived using bulk small RNA-seq from mouse cerebral cortices at E14, E17 and P0. a Bar plot showing the number of miRNAs in each co-expression network module. b Pearson correlation between developmental stages of the mouse cerebral cortex (E14, E17 and P0) and the module eigengenes (ME). ***p < 0.001, **p < 0.01, *p < 0.05. c, d Module eigengene values at each developmental stage for the black and green module. n = 10 biological replicates at E14, n = 12 in all remaining groups. e, f Top 15 enriched GO terms for the biological process (BP) category for the gene targets of the black and green modules.To identify key miRNAs at early and late embryonic brain developmental stages, we reconstructed the network structure of miRNAs within the black and green modules and identified key driver (hub) miRNAs as described previously34,35. The black module contained 21 hub miRNAs. 12 miRNAs were also upregulated in NPCs compared to neurons including members of the miR-92 family – miR-92a, miR-92b (Fig. 4c). Moreover, we detected 18 key driver miRNAs in the green module. Ten of these were also significantly upregulated in neurons compared to NPCs, including the neuron-specific miR-124 and miR-137 (Fig. 4d). This analysis confirmed the stage-specific regulation of a core set of miRNAs during embryonic cortical development.Fig. 4: Hub miRNAs in the black and green WGCNA modules.Network plot showing the hub miRNAs in the black (a) and green (b) modules. Hub miRNAs are represented as bigger nodes. Purple halos around the nodes indicate miRNAs that are significantly differentially expressed between E14 and P0. Differential expression of hub miRNAs from the black (c) and green (d) module detected in the NPCs versus neurons analysis. Positive log2 fold changes indicate an upregulation in NPCs, negative values correspond to upregulation in neurons. White bars indicate non-significant fold changes. e Distribution of differentially expressed genes in the mouse cerebral cortex between E14.5 and P0. The overall number of genes as well as the number of these genes that are validated targets from miRNAs in the black and green module as obtained from miRTarBase are given. The comparison of the distribution of differentially expressed targets of the black or green module from the overall distribution of differentially expressed genes was performed with a Fisher’s exact test. f Z-score transformed average expression of selected validated targets of miRNAs from the green module that are significantly upregulated at E14.5 compared to P0. The expression values in (e) and (f) were obtained from the bulk RNA-seq of the mouse cerebral cortex across development study by Weyn-Vanhentenryck et al.37.Importantly, all hub miRNAs in the green module and all but one miRNA in the black module were significantly differentially regulated at E14 compared to the P0 developmental stage (Fig. 4a, b). Furthermore, hub miRNAs in both the black and the green modules showed significantly stronger module membership and gene trait significance estimates compared to the remaining miRNAs in the respective module (Supplementary Fig. 3e–h), thereby confirming their importance. Here, module membership corresponds to the correlation between the miRNA’s expression and the module eigengene. Furthermore, the gene trait significance represents the correlation between miRNA expression and the trait of interest, in our study—the developmental stage36.To gain an insight into the relevance of hub miRNAs in regulating target gene expression across development, we re-analyzed bulk RNA-seq data from the mouse cortex at E14.5, E16.5, and P0 from the study by Weyn-Vanhentenryck and colleagues37. We limited this analysis to the experimentally validated targets of the conserved hub miRNAs that we obtained from miRTarBase38. We did not observe a significant deviation in the differential expression of target genes of the black module when comparing P0 and E14.5 (Fig. 4e). This could likely be explained by the overall lower number of validated targets for conserved key driver miRNAs from the black module (Supplementary Fig. 3d). However, the number of downregulated targets of the green module’s hub miRNAs was significantly higher compared to the overall background distribution of differentially expressed genes (Fig. 4e). The increased proportion of genes with a reduced expression at the time point when the green module’s hub miRNAs that target them have their expression peak strongly suggests that these miRNAs are involved in repressing these targets’ expression in vivo.To gain insight into the biological processes that the downregulated targets of green module hub miRNAs are involved in, we performed a GO annotation analysis (Fig. 4f). Notably, many of these genes are associated with the cell cycle, stem cell maintenance and proliferation (e.g., the transcription factors Sox2, Sox9 and Sox11), cell fate commitment (e.g., Neurog2), and negative regulation of neuronal differentiation. Furthermore, several targets are involved in activating the Notch signaling pathway (e.g., Notch1, Notch3) and negatively regulating the Wnt pathway. Important mediators of RNA splicing (e.g., the RNA binding proteins Ptbp1 and Ptbp2) were also among the downregulated targets. Notably, the green module hub miR-124-3p is known to repress Ptbp1 expression and thereby induce neuron-specific alternative splicing programs required for neuronal differentiation39. Therefore, these results indicate that miRNAs might be particularly important for the switch from undifferentiated NPCs to differentiated neurons by silencing the expression of key factors that promote stem cell maintenance, proliferation and non-neuronal splicing patterns.Prediction of miRNA co-targeting networks in the developing cerebral cortexWGCNA revealed that clusters containing multiple miRNAs are co-expressed during cortical development. Furthermore, when looking at the predicted targets of the black and green modules, we observed that only ~30% were potentially targeted merely by a single miRNA. On average, each gene was associated with more than 3 miRNAs (Supplementary Fig. 3a, b). This apparent redundancy of simultaneously expressed miRNAs that putatively bind to overlapping sets of genes supports the previously proposed hypothesis that miRNAs might act together to co-operatively exert stronger gene silencing (Fig. 1). To construct a comprehensive co-targeting network of miRNAs in the developing cerebral cortex, we designed a statistical framework that allowed us to detect which miRNAs share significantly more targets than expected by chance. In this analysis, we included all miRNAs whose module eigengene was correlated significantly with developmental time in the WGCNA (Fig. 3b). Furthermore, we only considered conserved miRNAs with at least 300 conserved targets and distinct seed sequences20. miRNAs belonging to the same broad conserved family and thereby having identical seed sequences, were grouped together for the statistical testing. After applying these filtering criteria, we performed the co-targeting prediction with a set of 77 miRNA families corresponding to 106 individual miRNAs (Fig. 5a). We detected 1216 significant co-targeting pairs with an adjusted p-value < 0.05 (Fig. 6, Supplementary data 4) and each miRNA had on average 41 co-targeting relationships (Fig. 5b). Interestingly, the number of significant co-targeting relationships was positively correlated with the miRNA’s average expression level (Fig. 5c). This finding implies that important miRNAs might need to be produced at higher amounts to facilitate their involvement in multiple cooperative gene silencing interactions in the co-targeting network. Importantly, hub miRNAs such as miR-92a/b from the black module as well as miR-124 and miR-137 from the green module were among the miRNAs with the highest expression level and highest number of co-targeting relationships (Fig. 5c) indicating that these miRNAs might be master regulators in both co-expression and multi-targeting networks. As expected, the length of the 3’ UTR was positively correlated with the number of significant miRNA pairs targeting a gene (Fig. 5d). Interestingly, Nova1, was among the top 10 genes involved in the highest number of co-targeting associations. NOVA1 is a neuron-specific splice factor that is crucial for neuronal viability40 and regulates the alternative splicing of genes involved in synapse formation41. In line with this, we detected a significant Nova1 expression increase in cortical samples at P0 compared to E14.5 (log2 fold change = 0.84, adj. p-value < 0.0001) when re-analyzing the RNA-seq data from Weyn-Vanhentenryck et al.37. Thus, our findings further highlight the potential importance of miRNAs in regulating neurodevelopment by influencing the expression of crucial effectors of alternative splicing processes.Fig. 5: Features of the miRNA co-targeting network of the developing cerebral cortex.The co-targeting network was constructed using miRNAs expressed in the mouse cerebral cortex at E14, E17 and P0. a Bar plot showing the number of miRNAs from each WGCNA module that were included in the co-targeting analysis after applying the filtering criteria. b Histogram of the number of miRNAs having a specific number of significant co-targeting relationships. c Scatter plot indicating the correlation between the number of co-targeting relationships with the log2-transformed miRNA’s average expression level; r corresponds to Pearson’s correlation coefficient. Labels indicate hub miRNAs from the WGCNA analysis of the black and green modules d Correlation between the 3’ UTR length of target genes and the number of significant co-targeting miRNA pairs the genes are targeted by; r corresponds to Pearson’s correlation coefficient. e Bar plot with the association of the number of intra-modular (between miRNAs of the same module) and inter-modular (involving miRNAs from different modules) co-targeting relationships. Odds ratios were obtained with a Fisher’s exact test. Error bars indicate the 95% confidence interval. f Heatmap with the top 10 strongest co-targeting relationships between pairs of miRNAs/conserved miRNA families. Values indicate odds ratios as measures of association between the co-targeting pairs.Fig. 6: Significant co-targeting miRNA pairs in the mouse cerebral cortex.Colored squares on the heatmap indicate significant co-targeting relationships. miRNAs that belong to the same conserved family and therefore have identical co-targeting relationships to each other, are indicated with gray rectangles on the vertical axis of the heatmap.Remarkably, we detected the strongest co-targeting relationship between miR-137 and the miR-25/miR-363/miR-92 family, sharing 353 common gene targets (Fig. 5f). While miR-137 was upregulated at later developmental stages, all members of the miR-25/miR-363/miR-92 cluster were regulated in the opposite direction. This prominent association between hub miRNAs, that were not co-expressed, prompted us to investigate the distribution of intra-modular (between miRNAs from the same module) and inter-modular (miRNAs from different modules) co-targeting relationships for the black and green WGCNA modules. While we expected that relationships between co-expressed miRNAs might be enriched, we indeed observed a comparable distribution of intra- and inter-modular co-targeting pairs (Fig. 5e). These findings point to two distinct potential regulatory mechanisms justifying the need for miRNAs to share significantly more targets than expected by chance. On the one hand, miRNAs with common targets regulated in opposite directions during cortical development might arise from the need to facilitate target gene regulation at different time points or in distinct cell types. On the other hand, co-expressed miRNAs that regulate the same genes could cooperatively bind to their targets to exert a stronger repressive effect compared to single miRNAs.Enhanced gene silencing effect by cooperative binding of co-expressed miRNAsTo validate the hypothesis that co-expressed miRNAs bind co-operatively to their common targets to enhance their repressive effect, we employed luciferase reporter assays focusing on miRNAs that were upregulated during embryonic brain development. We performed a gene ontology term look-up of all targets and filtered potential candidate genes for terms including brain/forebrain development and (central) nervous system development. Then we looked for genes with two or more conserved binding sites in their 3’ UTR for miRNAs upregulated during cortical development. Using this strategy, we compiled a list of 7 candidate target genes (Neurod1, Apc, Dcx, Ndst1, Zeb2, Src and Cxcl12) and 5 miRNAs (miR-128-3p, miR-129-5p, miR-135b-5p, miR-137-3p and miR-153-3p).First, we performed RT-qPCR analyses to validate the expression patterns that we observed in the bulk sequencing of the miRNAs selected for luciferase assays. All 5 miRNAs showed significantly increased expression levels in neurons compared to NPCs with mean relative quantification (RQ) values ranging from 2.9 to 3987 in the case of miR-135b-5p (Supplementary Fig. 4a–e). Furthermore, all miRNAs were associated with increased expression at later developmental stages in cortical tissue samples. The relative expression of miR-128-3p and miR-135b-5p significantly increased over all time points analyzed (Supplementary Fig. 4f, h). miR-137-3p and miR-153-3p were associated with a significantly higher expression levels at E17 and P0 compared to E14 (Supplementary Fig. 4 i, j). miR-129-5p was the only miRNA that showed a non-significant trend of higher expression at later developmental stages (Supplementary Fig. 4g).Subsequently, for each luciferase reporter assay, we cloned the target gene’s 3’UTR fragment containing the binding sites of the selected miRNAs, downstream of the synthetic Renilla luciferase gene into the psiCHECKTM-2 vector, which also expresses firefly luciferase as internal control. Plasmids were then co-transfected into HEK293 cells with different combinations of miRNA mimics and luciferase activity was evaluated in cell lysates 48 h later. The fold change of Renilla/firefly luciferase signal was normalized to a control group containing either a non-targeting siRNA or a miRNA mimic without a predicted binding site in the 3’UTR fragment of the target gene.The first gene we investigated, Neurod1, is an important transcription factor regulating neuronal differentiation and migration in the developing cerebral cortex42. Neurod1 contains binding sites for miR-137-3p and miR-153-3p that were previously demonstrated to co-operatively repress the gene’s expression20. Therefore, we used this experiment as a positive control. We independently cloned two Neurod1 3’ UTR fragments with different lengths (274 bp and 1188 bp, respectively). The luciferase activity was significantly reduced in cells transfected with the combination of miR-137-3p and miR-153-3p compared to individual miRNAs as well as two different negative controls (siRNA and miR-128-3p), thereby confirming the co-targeting effect previously observed (Fig. 7a, b).Fig. 7: Luciferase activity in lysates of HEK293 cells transfected with plasmids containing 3’ UTR fragments of target genes.a–d Lysates were co-transfected with different combinations of miRNA mimics. Fold change (FC) of luciferase activity was obtained by calculating the ratio of Renilla luciferase and firefly luciferase activity and then normalizing to the mean of the control siRNA or miRNA group. Locations of the binding sites in the 3’ UTR of the respective gene are represented by colored dots. The length of the 3’ UTRs is indicated in kbp. e–h Fold change of luciferase activity was measured after introducing point mutations in the binding sites of the miRNAs in the 3’ UTR regions. Colored dots on the x-axis indicate that the binding site of the respective miRNA was intact, gray dots indicate a mutated binding site. i–k Expression of miRNAs and target genes during embryonic brain development. Expression values were normalized to a 0–1 range. miRNA expression was quantified by bulk RNA-seq of the mouse cerebral cortex at E14, E17 P0. Target gene expression values were measured at E14.5, E16.5 and P0 and were obtained from the bulk RNA-seq study of the mouse cerebral cortex at different developmental stages by Weyn-Vanhentenryck et al.37. Data are shown as mean ± standard error of the mean and individual values. *p < 0.05, **p < 0.01, ***p < 0.001, two-way analysis of variance followed by Tukey’s post-hoc test. N = 12 replicates in (a) and (c), n = 18 in d, n = 9 in (b), (c) and (e–h), n = 2 in (i–k).Src, a proto-oncogene that codes for a non-receptor protein tyrosine kinase is another important neurodevelopmental regulator that contains binding sites for both miR-137-3p and miR-153-3p in its 3’UTR. Src is expressed at steadily high levels in the mouse neocortex from E12.5 to P1 and overexpression of this gene leads to impaired neuronal migration due to altered adhesion properties and cytoskeletal dynamics43. We observed that the combination of both miRNAs exerted a significantly stronger repressive effect on Src compared to each individual miRNA (Fig. 7c).As an additional target gene for the reporter assays, we selected Cxcl12 which is also an important regulator of early brain development. Cxcl12 is involved in migration of NPCs, early localization of Cajal-Retzius cells in the developing cortex, and in axon guidance and pathfinding44,45. Cxcl12 is a predicted target of miR-137-3p and miR-135b-5p. Luciferase reporter assays showed that co-transfection with both miRNAs significantly increased the repressive effect compared to the negative siRNA control and treatment with each miRNA separately (Fig. 7d).In contrast, we did not observe a significant co-targeting effect for miR-153-3p with miR-129-5p and miR-137-3p with miR-128-3p in reporter assays with Apc and Ndst1, respectively (Supplementary Fig. 5a, c). Furthermore, we aimed to investigate if higher-order co-targeting interactions, including three miRNAs would exert stronger silencing effects on gene expression. To this end, we co-transfected cells with plasmids containing the 3’ UTR fragment of the Dcx gene and mimics of miR-128-3p, miR-129-5p, and miR-135b-5p (Supplementary Fig. 5b). In an additional assay, we cloned the 3’ UTR of Zeb2 and co-transfected cells with its targeting miRNAs – miR-129-5p, miR-137-3p and miR-153-3p (Supplementary Fig. 5d). Treatment with different miRNA mimics significantly reduced luciferase activity compared to the siRNA control in the Zeb2 assays. However, we did not observe a significant cooperative repressive effect for either Dcx or Zeb2.To confirm that the cooperative silencing effect for Neurod1, Src, and Cxcl12 was mediated via direct binding of the miRNAs, we employed in vitro site-directed mutagenesis to introduce mutations in the miRNA binding sites of the 3’ UTRs. We generated plasmids carrying a mutation in one of the binding sites or in both binding sites of the respective co-targeting pair. For all three genes, we observed the highest luciferase activity in the cell lysates where both miRNA binding sites were mutated. Conversely, the strongest repressive effect was detected when both binding sites were intact, thereby confirming that miRNAs acted co-operatively to exert stronger gene silencing (Fig. 7e–h). Furthermore, we investigated the expression patterns of these three genes by employing the data from Weyn-Vanhentenryck et al.37. Neurod1 was significantly downregulated at P0 compared to E14.5 and E16.5, correlating negatively with the expression pattern of the co-targeting pair miR-137-3p and miR-153-3p (Fig. 7i, Supplementary Fig. 6a). While Src was not differentially regulated, it still had its expression minimum at P0 (Fig. 7j, Supplementary Fig. 6b). Interestingly, Cxcl12 had a significantly reduced expression only at E16.5 compared to E14.5, but this pattern negatively correlated with the expression peak of miR-135b-5p at E17 (Fig. 7k, Supplementary Fig. 6c). These matching in vivo expression patterns of the miRNAs and the genes they co-target lend further support to the biological relevance of the co-operative miRNA gene silencing we observed in the luciferase assays.
