A multifunctional locus controls motor neuron differentiation through short and long noncoding RNAs

Abstract

The transition from dividing progenitors to postmitotic motor neurons (MNs) is orchestrated by a series of events, which are mainly studied at the transcriptional level by analyzing the activity of specific programming transcription factors. Here, we identify a post-transcriptional role of a MN-specific transcriptional unit (MN2) harboring a lncRNA (lncMN2-203) and two miRNAs (miR-325-3p and miR-384-5p) in this transition. Through the use of in vitro mESC differentiation and single-cell sequencing of CRISPR/Cas9 mutants, we demonstrate that lncMN2-203 affects MN differentiation by sponging miR-466i-5p and upregulating its targets, including several factors involved in neuronal differentiation and function. In parallel, miR-325-3p and miR-384-5p, co-transcribed with lncMN2-203, act by repressing proliferation-related factors. These findings indicate the functional relevance of the MN2 locus and exemplify additional layers of specificity regulation in MN differentiation.

Keywords: competing endogenous RNA; long noncoding RNAs; microRNAs; motor neuron differentiation; single-cell RNA sequencing.

Figure EV1. LncMN2: locus structure, KO mutant generation, and analysis. Relative to Fig1

1. Directed MN differentiation of mESCs carrying a Hb9::GFP transgene. Timeline showing the days of treatment (green blocks) and the differentiation stages. Time‐points relevant in this study are highlighted: mESC (day0), MN‐committed EB (day6, EB6), and mature neural population cultured for 14 days after EB6 dissociation (DIV14).

2. qRT‐PCR analysis of cell type‐specific markers along MN differentiation: Nanog (mESC), Pax6 (MN progenitors), Hb9, and ChAT (postmitotic MNs). Differentiation stages (mESC and EB at specific maturation days) are indicated below. Expression peaks are set as 1. Expression levels (fold relative to the peak) are normalized over Atp5o mRNA.

3. 5330434G04Rik locus (here renamed as lncMN2) as visualized by the Ensembl genome browser. The lncRNA splicing isoforms and the embedded miR‐325‐3p and miR‐384‐5p are shown.

4. qRT‐PCR analysis of lncMN2‐202/203/204 isoforms (left), miR‐325‐3p and miR‐384‐5p (right) along MN differentiation. Data are normalized over Atp5o mRNA levels (for lncRNA isoform quantification) or over snoRNA U25 (for miRNA quantification).

5. qRT‐PCR analysis of lncMN2‐203 and lncMN2‐202/204 isoform in nuclear and cytoplasmic RNA extracts from EB6. Gapdh and snoRNA U25 are used as fractionation controls. Normalizations were performed on the total amount of RNA. Data, which are expressed as mean (error bars represent SD), represent the subcellular fraction of each target RNA over its total expression level. N = 3 biological replicates., *P ≤ 0.05, **P ≤ 0.01, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test).

6. Block‐and‐line model of the genome‐editing strategy to inhibit lncMN2 transcription. Schematic representation of the donor vector is shown. Homology arms (HAL and HAR) poly(A)/2×MAZ cassette (PAS) selection cassette (RFP‐T2A‐Puro‐PolyA) and LoxP sites are indicated. Cells were transfected with Cre recombinase to remove the selection cassette.

Figure 1. LncMN2 KO affects MN differentiation

1. Block‐and‐line model of lncMN2 gene structure. Exons specific for lncMN2‐202/204 isoforms are in yellow; lncMN2‐203‐specific exon is in blue. Shared exons are gray. miR‐325 and miR‐384 are represented as red and green blocks, respectively.

2. qRT‐PCR expression analysis of lncMN2 isoforms (left panel) and embedded miRNAs (right panel) in mESC and EBs differentiated for 6 days (EB6). Data, that are expressed as mean (error bars represent SD), were normalized over Atp50 mRNA (left) and snoRNA U25 (right) levels and reported as RNA fold over the control (mESC), set as 1. N = 3 biological replicates.

3. qRT‐PCR expression analysis of lncMN2‐202/203/204 isoforms (left panel), miR‐325‐3p (middle panel), and miR‐384‐5p (right panel) in WT, KO1, and KO2 EB6. Details as in B.

4. Cytofluorimetric analysis of cells from dissociated EB6. The left diagram reports the distribution of cells (percentage of the peak) according to GFP levels, from WT (gray), KO1 (orange), and KO2 (blue) EB. In the histogram aside, the percentage of GFP⁺ cells in WT, KO1, and KO2 is reported. Data are expressed as mean and error bars represent SD. N = 3 biological replicates.

5. Left: representative immunofluorescence of the MN marker Isl1 (red signals) in cells maintained in vitro for 14 days after EB6 dissociation (DIV14). DAPI counterstaining in blue. Scale bar = 100 µm. Right: box plot shows the distribution of Isl1⁺ nuclei with respect to the total cell number marked with DAPI staining in WT and in lncMN2 KO mutants (clone 1 and clone 2). Data are expressed as mean and error bars represent SD. The central band represents the median and the whiskers represent the distribution that goes from the minimum value to the lower quartile and from the upper quartile to the maximum value. N = 3 biological replicates and 10 fields were analyzed for each replicate.

Data information: Histogram dots represent single fields replicates. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test). Source data are available online for this figure.

Figure EV2. LncMN2: generation of ∆miRNAs mutants. Relative to Fig1

1. Block‐and‐line model of miR‐325 deletion in lncMN2 gene. The deleted region is indicated as ∆miR‐325. Further details as in Fig 1A.

2. Alignments of lncMN2‐203 WT (Seq_2) and ∆325 (Seq_1) sequences highlight a 124‐bp‐long genomic deletion.

3. Block‐and‐line model of miR‐384 deletion in lncMN2 gene. The deleted region is indicated as ∆miR‐384. Further details as in Fig 1A.

4. Alignments of lncMN2‐203 WT (Seq_2) and ∆384 (Seq_1) sequences highlight a 215‐bp‐long genomic deletion.

5. Left panel. qRT‐PCR expression analysis of miR‐325‐3p in ∆325 mESCs (clones 1 and 2) differentiated for 6 days (EB6). Data, which are expressed as mean (error bars represent SD), are normalized over snoRNA U25 levels and reported as RNA fold over the control (WT EBs), set as 1. Right panel. Percentage of GFP+ cells (MNs) in WT and ∆325 (clones 1 and 2) upon cytofluorimetric analysis of cells from dissociated EB6. N = 3 biological replicates.

6. Left panel. qRT‐PCR expression analysis of miR‐384‐5p in ∆384 mESCs (clones 1 and 2) differentiated for 6 days (EB6). Data, which are expressed as mean (error bars represent SD), are normalized over snoRNA U25 levels and reported as RNA fold over the control (WT EBs), set as 1. Right panel. Percentage of GFP+ cells (MNs) in WT and ∆384 (clones 1 and 2) upon cytofluorimetric analysis of cells from dissociated EB6. N = 3 biological replicates.

Data information: **P ≤ 0.01, ***P ≤ 0.001, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test).

Figure 2. miR‐466i‐5p interacts with lncMN2‐203

1. A

   Block‐and‐line model of the lncMN2‐203 structure. Details as in Fig 1A. The miRNA‐interacting region in exon5 (dark blue) is represented as a light blue block. Sequences of lncMN2‐203 predicted miRNA‐binding sites and of miR‐466i‐5p (above) or miR‐669a‐5p (below) are indicated in the textbox.

2. B

   qRT‐PCR expression analysis of miR‐466i‐5p and miR‐669a‐5p in mESC and EB6. Data, which are expressed as mean (error bars represent SD), are normalized over snoRNA U25 levels and reported as RNA fold over the control (mESC), set as 1. N = 3 biological replicates.

3. C

   Schematic representation of the luciferase‐based reporter constructs. 1,500 bps of WT lncMN2‐203 exon5 (dark blue) was cloned downstream to the Renilla luciferase ORF, represented in gray (Luc‐WT). The predicted miRNA‐binding region is light blue. In the ∆SP mutant (Luc‐∆SP), a 260 bp deletion was introduced. As indicated, each construct was co‐transfected with LNAs targeting miR‐466i‐5p or miR‐669a‐5p.

4. D

   Quantification of Renilla luciferase activity in NSC‐34 cells co‐transfected with the Luc‐WT or the Luc‐∆SP construct and LNA against miR‐466i‐5p or miR‐669a‐5p. The histogram represents the mean luciferase activities (error bars represent SD) folded over the control, set as 1 (basal activity from the Luc‐WT in the presence of scrambled LNA). N = 3 biological replicates.

5. E

   Ago2 cross‐linked RNA immunoprecipitation assay performed on total extracts from WT EB6. qRT‐PCR quantification of lncMN2 enrichments (isoform −203 or −202/204) in the IP and IgG samples, relative to the input is shown. Ptbp1 and Gapdh mRNAs were used as positive and negative controls, respectively. Values are expressed as input percentages and error bars represent SD. N = 3 technical replicates.

6. F–I

   RNA pull‐down assay performed on total extracts from differentiated WT EB6. qRT‐PCR quantification of lncMN2‐203 (F), miR‐466i‐5p (G), miR‐669a‐5p (H) or miR‐325‐3p (I) enrichments in the lncMN2‐203 specific (MN2) or in the control (U1) pull‐down samples. Values are expressed as input percentages and error bars represent SEM. N = 3 biological replicates.

Data information: Histogram dots represent single replicates. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test).

Figure EV3. LncMN2‐203 potentially binds miRNAs. Relative to Fig 2

1. Left: LncMN2‐203 secondary structure as predicted by the Mfold. Right: Magnification of the miRNA‐binding region.

2. Western blot analysis of Ago2 immunoprecipitation from CLIP assay performed on EB6 total extracts. Immunodetection of Ago2 in Input (Inp) and immunoprecipitated (IP) protein fractions is shown.

3. RNA pull‐down assay performed on total extracts from differentiated WT EB6. qRT‐PCR quantification of U1 enrichment in the lncMN2‐203‐specific (MN2) or in the control (U1) pull‐down samples. Values are expressed as input percentages, and error bars represent SEM. N = 3 biological replicates.

Figure 3. LncMN2 ∆SP mutation impairs MN production

1. Block‐and‐line model of the lncMN2 edited locus. Details as in Figs 1A and 2A. Predicted miRNA‐binding site deletion is represented as ∆SP.

2. qRT‐PCR expression analysis of lncMN2‐203 isoform, lncMN2‐202/204, miR‐325‐3p, and miR‐384‐5p embedded miRNAs in WT, ∆SP1, and ∆SP2 EB6. Data, which are expressed as mean (error bars represent SD), are normalized over Atp5o mRNA levels (for lncRNA isoform quantification) or over snoRNA U25 levels (for miRNA quantification) and are reported as RNA fold over the control (mESC), set as 1. N = 3 biological replicates.

3. Cytofluorimetric analysis of cells from dissociated EB6. The right diagram reports the cell distribution (percentage of the peak) according to GFP levels, from WT (gray line), ∆SP1 (purple line), and ∆SP2 (blue line) EB6. In the histogram aside, the percentage of GFP⁺ cells in WT, ∆SP1, and ∆SP2 is reported. Data are expressed as mean, and error bars represent SD. N = 3 biological replicates.

4. Left: representative immunofluorescence of the MN marker Isl1 (red signals) in cells maintained in vitro for 14 days after EB6 dissociation (DIV14). DAPI counterstaining in blue. Scale bar = 100 µm. Right: box plot shows the distribution of Isl1⁺ nuclei with respect to the total cell number marked with DAPI staining in WT and in ∆SP mutants (clone 1 and clone 2). The central band represents the median, and the whiskers represent the distribution that goes from the minimum value to the lower quartile than from the upper quartile to the maximum value. N = 3 biological replicates and 10 fields were analyzed for each replicate.

Data information: Histogram dots represent single‐field replicates. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test). Source data are available online for this figure.

Figure EV4. Cell subpopulation marker expression and trajectory. Relative to Fig 4

1. Left: UMAP plot of the integrated dataset depicting cell identity assignment to NP (violet dots), MNP (green dots), EMN (orange dots), LMN (red dots), IN (light blue), or NA (gray dots) subpopulations; Right: UMAP plot of the integrated dataset depicting cell identity assignment to proliferating (blue dots) or postmitotic (red dots) subpopulations.

2. Single‐cell expression of marker genes over the UMAP representation of the integrated dataset.

3. UMAP and related cell density plot of the whole dataset with Monocle software trajectory. The WT, ∆SP, and KO samples are shown in red, green, and blue, respectively.

4. Dot plot generated by the Seurat DotPlot function represents the expression of cell identity markers (listed on the left) for each identified subpopulations in every condition (WT, ∆SP and KO). Genes used for the assignment of cluster identity are marked in red. Dataset conditions and cell types are indicated on the bottom.

Figure 4. SC‐seq reveals the implication of lncMN2 in motoneurogenesis

1. UMAP visualization plots for WT2 (1,435 cells), ∆SP1 (2,039 cells), and KO1 (1,955 cells) produced by Monocle3 (Trapnell et al, 2014). The plots present the estimated trajectories (dark gray), and the colors correspond to the cell identity: neural precursors NP cells (purple), motor neurons progenitors (MNP, green), early motor neurons (EMN, orange), late motor neurons (LMN, red), and interneurons (IN, light blue). The cells not assigned to any of the previous groups (NA) are indicated in gray.

2. Histograms reporting the percentage of cells detected as proliferating and postmitotic (Dataset EV1). All percentages of cells were calculated, for each condition, on samples belonging to the same batch (WT2, ∆SP1, ∆SP2, KO1, and KO2) excluding not assigned cells. The significance of cell subpopulation differences between WT and mutants was calculated with the Fisher's exact test (see Materials and Methods). *P < 0.05, **P < 0.01, ***P < 0.001, n.s. > 0.05; N = 2 biological replicates.

3. Histograms reporting the percentage of cells assigned to NP, MNP, EMN, LMN, and IN (Dataset EV1). All percentages of cells were calculated, for each condition, on samples belonging to the same batch (WT2, ∆SP1, ∆SP2, KO1, and KO2) excluding not assigned cells. The significance of cell subpopulation differences between WT and mutants was calculated with the Fisher's exact test (see Materials and Methods). *P < 0.05, **P < 0.01, ***P < 0.001, n.s. > 0.05. N = 2 biological replicates.

4. Heatmap generated by the Seurat DoHeatmap function (Stuart et al, 2019) represents the expression of cell identity markers (listed on the left) for each cell of the identified subpopulations. Genes used for the assignment of cluster identity are marked in red. Dataset samples are indicated on the top. Red represents the maximum expression value, light blue the minimum.

Figure 5. GO enrichment analyses in MNP and LMN

1. Volcano plots representing differentially expressed genes between WT MNP and ∆SP MNP (left) and WT MNP and KO MNP (right). Each dot represents a gene; −log₁₀ P value is represented in y‐axis while average log₂FC in x‐axis. Significantly deregulated genes are depicted as black dots.

2. Volcano plots representing differentially expressed genes between WT LMN and ∆SP LMN (left) and WT LMN and KO LMN (right). Each dot represents a gene. −log₁₀ P‐value is represented in y‐axis while average log₂FC in x‐axis. Significantly deregulated genes are depicted as black dots.

3. GO‐enriched categories in MNP related to downregulated (left) and upregulated (right) genes in ∆SP (top) and KO (bottom) mutants. Green dots indicate “cell cycle” category.

4. Bar plot representing the number of “cell cycle”‐related genes upregulated in both ∆SP and KO MNP (light green) or specifically in either one of the two mutants (dark green).

5. GO‐enriched categories in LMN cells related to downregulated (left) and upregulated (right) genes in ∆SP (top) and KO (bottom) mutants. Red dots indicate “nervous system development” category.

6. Bar plot representing the number of “nervous system development”‐related genes downregulated in both ∆SP and KO (dark red) or specifically in either one of the two mutants (light red).

Figure EV5. GO enrichment analyses in EMN. Relative to Fig 5

1. Volcano plots representing differentially expressed genes between WT EMN and ∆SP EMN (left) and WT EMN and KO EMN (right). Each dot represents a gene. ‐log₁₀ P‐value is represented in y‐axis while average log₂FC in x‐axis. Significantly deregulated genes are depicted as black dots.

2. GO‐enriched categories in EMN related to downregulated (left) and upregulated (right) genes in ∆SP (top) and KO (bottom) mutants. Green dots indicate “cell cycle” category.

Figure 6. Neuronal differentiation genes respond to mir‐466i‐5p activity

1. Venn diagram showing the relationships between genes upregulated in KO and ΔSP mutants. Gene numbers and percentages of overlap are shown below diagrams for both groups. The intersection between miR‐325‐3p and miR‐384‐5p predicted targets and genes specifically upregulated in KO mutants are indicated in green (MNP, left) or in red (LMN, right).

2. qRT‐PCR expression analysis of Fbxw11, H2afz, and Phb2 in NSC‐34 transfected with miR‐325‐3p LNA and Sapcd2, Suv39h, and Tdrkh in NSC‐34 transfected with miR‐384‐5p LNA. Data, which are expressed as mean (error bars represent SD), are normalized over Atp5o mRNA levels and reported as RNA fold over the control (scrambled LNA, CTRL), set as 1. N = 3 biological replicates.

3. Venn diagram showing the relationships between genes downregulated in KO and ∆SP mutants. Gene numbers and percentages of overlap are shown below diagrams for both groups. The intersection between miR‐466i‐5p‐predicted targets and genes downregulated in both mutants are indicated in the MNP (left) and LMN (right) subpopulations as the green and red area, respectively.

4. Bar plot depicting for each analyzed cell subpopulation (MNP, EMN, and LMN) the percentage of miR‐466i‐5p targets among genes downregulated in both ∆SP and KO mutants. Statistical significance was assessed using the Fisher's exact test. *P ≤ 0.05, **P ≤ 0.01, n.s. > 0.05.

5. Bar plot depicting the ratio of miR‐466i‐5p predicted targets (x‐axis) in GO‐enriched categories of genes downregulated in both ∆SP LMN and KO LMN. Neural categories (“modulation of chemical synaptic transmission” and "nervous system development”) are colored in red. GO enrichment analysis was performed using WebGestalt software.

6. Violin plots representing for each condition (WT, ∆SP, and KO) the expression of Scrt1, Ncan, Nrxn1, and Snap25 in LMN. Plots were produced by Seurat “VlnPlot” function. N = 2 biological replicates. Statistical significance was assessed using the Wilcoxon rank‐sum test (default Seurat FindMarkers function). ***P ≤ 0.001.

7. qRT‐PCR expression analysis of Scrt, Ncan (up) and Nrxn1, Snap25 (down) in NSC‐34 cells transfected with miR‐466i‐5p LNA. Data, which are expressed as mean (error bars represent SD), are normalized over Atp5o mRNA levels and reported as RNA fold over the control (scrambled LNA, CTRL), set as 1. N = 3 biological replicates.

8. qRT‐PCR expression analysis of Scrt, Ncan (up) and Nrxn1, Snap25 (down) in ∆SP mESCs transfected with miR‐466i‐5p LNA and differentiated to early EBs (day 3). Data, which are expressed as mean (error bars represent SD), are normalized over Atp5o mRNA levels and reported as RNA fold over the control (scrambled LNA, CTRL), set as 1. N = 3 biological replicates.

Data information: Histogram dots represent single replicates. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, n.s. > 0.05 (two‐tailed, unpaired Student’s t‐test).

Figure 7

Schematic representation of the circuitry linking lncMN2‐203, miR‐466i‐5p, mir‐325‐3p, miR‐384‐5p, and MN differentiation. The model was created with BioRender.com