Identifying key underlying regulatory networks and predicting targets of orphan C/D box SNORD116 snoRNAs in Prader-Willi syndrome

Abstract

Prader-Willi syndrome (PWS) is a rare neurodevelopmental disorder characterized by neonatal hypotonia, followed by hyperphagia and obesity. Most PWS cases exhibit megabase-scale deletions of paternally imprinted 15q11-q13 locus. However, several PWS patients have been identified harboring much smaller deletions encompassing the SNORD116 gene cluster, suggesting these genes are direct drivers of PWS phenotypes. This cluster contains 30 copies of individual SNORD116 C/D box small nucleolar RNAs (snoRNAs). Many C/D box snoRNAs have been shown to guide chemical modifications of RNA molecules, often ribosomal RNA (rRNA). Conversely, SNORD116 snoRNAs show no significant complementarity to rRNA and their targets are unknown. Since many reported PWS cases lack their expression, it is crucial to identify the targets and functions of SNORD116. To address this we modeled PWS in two distinct human embryonic stem cell (hESC) lines with two different sized deletions, differentiated each into neurons, and compared differential gene expression. This analysis identified a novel set of 42 consistently dysregulated genes. These genes were significantly enriched for predicted SNORD116 targeting and we demonstrated impacts on FGF13 protein levels. Our results demonstrate the need for isogenic background comparisons and indicate a novel gene regulatory network controlled by SNORD116 is likely perturbed in PWS patients.

Graphical Abstract

Figure 1.

Summary of genomic editing and characterization of differentiated neurons. (A) A schematic of the deletions present in our model system. Not drawn to scale. (B) UCSC Browser image of the chromosome 15q11-q13 locus displaying representative bigwig tracks from each genetic background and genotype. Blue tracks show RNA signal from the sense strand; red tracks show RNA signal from the antisense strand. Top track shows CRISPR gRNA binding locations; gray shading indicates deleted region. GENCODEv25 gene annotations are shown at the bottom; protein-coding genes are shown in blue, noncoding genes are shown in green, and To Be Experimentally confirmed (TEC) biotype genes are shown in red. Some isoforms have been removed for clarity. (C) Box and whisker plots displaying a subset of significant DEGs (p.adjust < 0.05) in all cell lines across both genetic backgrounds as inducible neurons compared to wild type H9 ESCs (n = 4–6 biological replicates). ESC samples are shown in black and neuron samples are shown in blue or red. Genotype is represented by data point shape. The y-axis represents log₂(foldChange) and the x-axis displays individual gene names for either pluripotency markers (left) or neuronal markers (right). Pseudocount was added to counts of all genes prior to calculation of log₂(foldChange). (D) Representative immunofluorescent imaging of CT2 lgDEL cell line for canonical neuronal proteins TUBB3 and MAP2 in addition to nuclear marker DAPI at day 11 post-induction. Images were acquired at 63x magnification. Scale bar represents 25 μm.

Figure 2.

Analysis of gene set perturbed in lgDEL cell lines versus controls as neurons. (A) Venn diagram displaying overlap of significant DEGs (p.adjust < 0.05) for lgDEL lines in both genetic backgrounds versus their isogenic WT controls. Left side of diagram represents significant downregulated (log₂FoldChange < 0) DEGs with shared genes highlighted in blue; right side of diagram represents significant upregulated (log₂FoldChange > 0) DEGs with shared genes highlighted in red. Significance of overlaps (P < 0.0001) determined via a permutation test. (B) Box and whisker plot showing expression of a subset of protein-coding genes in the chromosome 15q11-q13 region. Pseudocount was added to counts of all genes prior to calculation of log₂(foldChange). Significant DEGs (p.adjust < 0.05) are shown in orange. (C) Heatmap showing 50 most dysregulated significant DEGs in lgDEL vs WT H9 genetic background. Top 25 up- and downregulated genes were determined by average log₂(foldChange) between CT2 and H9 backgrounds. Shading indicates row z-score, with blue denoting downregulated gene expression and red denoting upregulated gene expression. Rows represent samples; columns represent individual genes. (D) Gene-concept network plot displaying GO terms of the molecular function (MF) category for all shared dysregulated genes. Main nodes (tan) correspond to the MF category with colored lines connecting to nodes of genes found in each category. Size of the main node corresponds to the number of DEGs in our data set contained within each ontology term. Colors of the gene nodes correspond to the log₂(fold Change) for each DEG in lgDEL vs WT control; red indicates log₂(foldChange) > 0, blue indicates log₂(foldChange) < 0. (E) Dot plot displaying disease ontology results for shared downregulated genes. The x-axis represents the log₂ fold enrichment value, and y-axis shows disease ontology terms. Size of the dot corresponds to the number of DEGs in our data set contained within each ontology term. Shading of the dot corresponds to the negative log₁₀ of the adjusted P-value.

Figure 3.

Analysis of gene set perturbed in smDEL versus controls as neurons. (A) Venn diagram displaying overlap of significant DEGs (p.adjust < 0.05) for smDEL lines in both genetic backgrounds versus their isogenic WT controls. Left side of diagram represents significant downregulated (log₂FoldChange < 0) DEGs with shared genes highlighted in blue; right side of diagram represents significant shared upregulated (log₂FoldChange > 0) DEGs with shared genes highlighted in red. Significance of overlaps (P < 0.0001) determined via a permutation test. (B) Box and whisker plot showing expression of a subset of protein-coding genes in the chromosome 15q11-q13 region. Pseudocount was added to counts of all genes prior to calculation of log₂(foldChange). Significant DEGs (p.adjust < 0.05) are shown in orange. (C) Heatmap showing 50 most dysregulated significant DEGs in smDEL versus WT H9 genetic background. Top 25 up- and downregulated genes were determined by average log₂(foldChange) between CT2 and H9 backgrounds. Shading indicates row z-score, with blue denoting downregulated gene expression and red denoting upregulated gene expression. Rows represent samples; columns represent individual genes. (D) Dot plot displaying disease ontology results for all shared dysregulated genes in smDEL lines across both backgrounds. The x-axis represents the log₂ fold enrichment value, and y-axis shows top 25 disease ontology terms. Size of the dot corresponds to the number of DEGs in our data set contained within each ontology term. Shading of the dot corresponds to the negative log₁₀ of the adjusted P-value. (E) Venn diagram displaying overlap of all significant DEGs (p.adjust < 0.05) after additional filtering for both genetic backgrounds and genotypes versus isogenic WT controls. Yellow shading denotes significant shared dysregulated gene set. Significance of overlap (p < 0.001) determined via a permutation test.

Figure 4.

Interrogation of shared dysregulated genes and analysis of SNORD116 predicted targets. (A) Dot plot displaying DisGeNET results for 42 shared dysregulated genes. The x-axis represents the log₂ fold enrichment value, and y-axis shows top 25 ontology terms. Size of the dot corresponds to the number of DEGs in our data set contained within each ontology term. Shading of the dot corresponds to the negative log₁₀ of the adjusted P-value. Note Tumor Cell Invasion, Liver carcinoma, Tumor Progression ontologies contain SNHG14 as one of the DEGs driving these categories. (B) Bar plot displaying the number of predicted targeting events per copy of SNORD116. Colors of the bars correspond to the three subgroups of SNORD116: SNORD116-I (copies 1–9), SNORD116-II (copies 10–24) and SNORD116-III (copies 25–30). (C) Plot displaying distribution of prediction interactions for SNORD116-III. The x-axis corresponds to the relative position within snoRNA copies, and y-axis represents the number of predicted interactions for which the center of the predicted binding interaction was used (black line). Color-coded bar on the x-axis indicates the position of C/C’ and D/D’ boxes found in snoRNA copies, indicated by green and purple, respectively. (D) Bar charts representing the proportion of exon, intron, and intron-exon junctions in the shared 42 dysregulated genes (Background: Shared Genes) and the predicted targeting of SNORD116-III copies on those shared genes (SNORD116-III vs Shared Genes). Exon category is subdivided based on genic location and displayed as donut plots. Coloring of donut plots is based on exon category; 5′UTRs are represented in orange, 3′UTRs are represented in blue, CDS is represented in yellow, and any portion of exonic sequence not falling under those categories is termed ‘other’ and shown in black. (E) Metagene plot of predicted binding sites for SNORD116 versus the shared dysregulated gene set. Black line shows average of all SNORD116 groups. Various pink lines show each individual group (I-III). Metagene coordinates (x-axis) of 0–1 represent the 5′UTR, coordinates 1–2 represent the gene body, and coordinates 2 + represent the 3′UTR. The density of the predicted binding is on the y-axis.

Figure 5.

Investigation of FGF13, a predicted SNORD116 target. (A) UCSC Browser image of the FGF13 locus displaying BEDtracks of SNORD116 predicted binding. Chromosome ideogram indicates location of FGF13 on X chromosome in red. Top track shows all FGF13 gene isoforms in GENCODEv25 annotation. Bottom track shows zoomed in view, with predicted binding shown to occur in 5′ exonic region of one transcript of FGF13 (ENST00000315930.10). (B) Cartoon depicting one predicted RNA-RNA interactions between SNORD116-25 and FGF13. Created in BioRender. Gilmore, R. (2024) BioRender.com/k76r157. Not drawn to scale. (C) Western blot image of FGF13 and GAPDH in CT2 WT and CT2 smDEL inducible neurons harvested at day 11 post-induction (n = 3 biological replicates each). Barplot shows quantification of FGF13 decrease. Significance (P = 0.05) determined by the Wilcoxon Rank Sum test.