Neurotypical control teeth were obtained through the Department of Pediatric Dentistry and Community Oral Health at the University of Tennessee Health Science Center (UTHSC). Teeth from children with PWS subtypes were collected remotely by the caregivers of these subjects after confirmation of the underlying genetic diagnosis. Subjects provided informed consent for tooth collection along with a Social Communication Questionnaire (SCQ) to assess ASD status. Tooth pulp was cultured from teeth and cell lines frozen during early passages in our DPSC Repository as previously described (Goorha and Reiter, 2017). The DPSC Repository and subsequent molecular studies were approved by the UTHSC institutional review board prior to conducting research.

Dental pulp stem cells (DPSC) used in this study were isolated and cultured according to our previously described protocol and stored in the DPSC Repository (Goorha and Reiter, 2017). Briefly, after mincing the dental pulp from inside the tooth cavity, 3 mg/mL Dispase II and 4 mg/mL Collagenase I were added to digest the tissue. Cells were then seeded on poly-D-Lysine coated 12-well plates with DMEM/F12 1:1, 10% fetal bovine serum (FBS), 10% newborn calf serum (NCS), and 100 U/mL penicillin and 100 ug/mL streptomycin (Pen/Strep) (Fisher Scientific, Waltham, MA). Confluent cultures (80%) were passaged with TrypLE^TM Express and neuronal differentiation performed only on early passage cells (<passage 4).

Dental pulp stem cells (DPSC) lines were seeded at 20,000 cells/cm² on poly-D-lysine coated plates or chamber slides (Ibidi, Planegg, Germany) with DMEM/F12 1:1, 10% fetal bovine serum (FBS), 10% newborn calf serum (NCS), with 100 U/mL penicillin and 100 ug/mL streptomycin (Pen/Strep). At 80% confluence, the neuronal differentiation protocol was followed as previously published in Kiraly et al. (2009) with an extended maturation phase (3 weeks versus 7 days) (Goorha and Reiter, 2017). Briefly, epigenetic reprogramming was performed by exposing the DPSC to 10 μM 5-azacytidine (Acros Scientific, Geel, Belgium) in DMEM/F12 containing 2.5% FBS and 10 ng/mL bFGF (Fisher Scientific, Waltham, MA) for 48 h. Neural differentiation was induced by exposing the cells to 250 μM IBMX, 50 μM forskolin, 200 nM TPA, 1mM db-cAMP (Santa Cruz, Dallas, TX), 10 ng/mL bFGF (Invitrogen, Carlsbad, CA), 10 ng/mL NGF (Invitrogen, Carlsbad, CA), 30 ng/mL NT-3 (Peprotech, Rocky Hill, NJ), and 1% insulin-transferrin-sodium selenite premix (ITS) (Fisher Scientific, Waltham, MA) in DMEM/F12 for 3 days. At the end of the neural induction period, neuronal maturation was performed by maintaining the cells in Neurobasal A media (Fisher Scientific, Waltham, MA) with 1mM db-cAMP, 2% B27, 1% N2 supplement, 30 ng/mL NT-3, and 1X Glutamax (Fisher Scientific, Waltham, MA) for 21 days.

Once DPSC-neurons were matured for 3 weeks, total RNA was collected using the Zymo Directzol RNA extraction kit (Zymo, Irvine, CA). Extracted RNA was assayed for quality and integrity using the Agilent Bioanalyzer 6000 pico chip (Agilent, Santa Clara, CA). Only RNA with an RNA Integrity Number (RIN) ≥ 9.0 was used for RNAseq studies. Library preparation and RNAseq was performed by Novogene (NovoSeq 6000) (Sacramento, CA) using the Illumina platform and paired end reads. 20M reads per sample were collected.

FASTQ files from Novogene were analyzed for quality and trimmed using FASTQC. All reads were trimmed to remove nucleotides with Phred scores < Q20. The trimmed FASTQ files were aligned to the human genome reference library hg19 using RNASTAR (Widmann et al., 2012). Once aligned, the SAM files were collected and mined for read count information of each gene present in the reference file. Read counts were normalized using Counts per Million (CPM) method across groups for the entire experiment. Principle component analysis and Pearson’s coefficient plots were performed on the normalized transcriptome profile. A Wilcoxon’s t-test was used to determine significance across groups. All genes that fail to yield a p-value ≤ 0.05 and a fold change greater than 1.5 were removed. Benjamini and Hochberg false discovery rate (FDR) was performed on this trimmed gene list. All genes that failed to yield an FDR rate of ≤ 0.05 were removed. The final significant differential gene lists were loaded into the ClustVis web tool to generate heatmaps (Metsalu and Vilo, 2015). Additionally, the targets were loaded into the web based enrichment analysis tool, Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang da et al., 2009b) to identify enriched gene ontology (GO) terms.

Dental pulp stem cell (DPSC) were grown and differentiated on 3-well chamber slides (Ibidi, Planegg, Germany) coated with poly-D-lysine. Cells were fixed using a 4% paraformaldehyde solution for 10 min. Once fixed, cells were blocked and permeabilized using PBS with 1% BSA, 10% FBS, and 0.3% Triton X-100 for 1 h. Primary antibodies were diluted to 1:500 for anti-Beta Tubulin (Millipore, ab9354) and anti-TOMM20 (Santa Cruz, sc-17764) in the blocking solution and incubated overnight at 2–8°C with agitation. After overnight incubation, the slides were washed 3x with PBS-T for 10 min before the secondary antibodies, Goat anti-mouse Alexa Fluor 488 (Life Technologies, A11029) and Goat anti-chicken Alexa Fluor 594 (Life Technologies, A11042) were added at a 1:1000 dilution. Slides were incubated at room temperature for 1 h and then washed again 3× for 10 min each. Finally, Prolong Gold Antifade with DAPI (Fisher Scientific, Waltham, MA) was applied for mounting. Slides were imaged on a Zeiss 710 confocal microscope at 63× magnification using Z-stacking to image the entire neuron via ZEN software (Black Edition).

Western blots were performed as previously described (Urraca et al., 2018). Briefly, protein was extracted from 3-week mature neurons using neuronal protein extraction reagent (N-PER) (Fisher Scientific, Waltham, MA) and protease inhibitor cocktail (Roche). Samples were resolved on a NuPage 1.5mm 4–12% Bis-Tris gel (Invitrogen, Carlsbad, CA) according to manufacturer instructions and transferred to an Invitrolon^TM -PVDF membrane (Invitrogen, Carlsbad, CA). The membrane was blocked using Odyssey Blocking Buffer (Licor, Lincoln, NE) for 1 h and incubated overnight at 2–8°C in primary antibody with agitation. Primary antibodies used: anti-MAP2 (Santa Cruz, sc-32791), anti-Nestin (Santa Cruz, sc-23927), and anti-GABA A receptor beta 3 (Abcam, ab98968). anti-GAPDH (Abcam, ab157156) was used as a protein loading control. Blots were incubated at room temperature for 1 h in secondary antibodies for both the 700 and 800 channels using Li-Cor IR secondary antibodies (Licor, 926-32212 and 926-68074). Blots were imaged on a Li-Cor Odyssey^TM Fc Imager. Both the 700 and 800 channels were exposed for 2 min.

All imaging analysis was done using coded cell lines, so the observer was blinded to the genotypes. Only after data collection and analysis were cell lines decoded. For each cell line (≥ 5 individuals per group), ≥15 neurons were imaged for analysis of mitochondrial abundance within the neuron. To quantify mitochondrial area, cells were labeled with anti-Beta Tubulin (total neuronal area) and anti-TOMM20 (mitochondrial area), images were loaded into the Imaris imaging analysis software (Oxford Instruments, Abingdon, United Kingdom) and surfaces for each marker were created (see Figure 6). The mitochondrial area (anti-TOMM20) was then divided by the total neuronal area (anti-Beta Tubulin) to give the percentage of the neuron that contains mitochondria. After collecting all data, significance testing was performed by both ordinary one-way ANOVA and Turkey’s multiple comparisons test with a single pooled variance.

Several studies have shown an incidence of about 35% for ASD in PW-UPD subjects compared to 18% in PW-del subjects (Bennett et al., 2015; Baker et al., 2018; Dimitropoulos et al., 2019). To measure ASD remotely, we used the Social Communication Questionnaire (SCQ). The SCQ provides a way to ascertain the likelihood a subject may have ASD, especially when face to face evaluation is not possible. Several groups have now validated the SCQ as a screening tool for ASD against other tools like the Autism Diagnostic Interview, Revised (ADI-R) (Hutchins et al., 2007; Sappok et al., 2015; Barnard-Brak et al., 2016; Chesnut et al., 2017; Mouti et al., 2019). We use the “lifetime” version of the test which has been proven to be a more accurate predictor of ASD (Chesnut et al., 2017). A cutoff score of 15 is considered “possible ASD” according the SCQ lifetime tool. We selected 4 lines from each group (neurotypical control, PW-del, PW-UPD -ASD, and PW-UPD + ASD) for RNAseq analysis and additional lines were used for the immunofluorescence analysis. Figure 1 shows the SCQ scores from these subjects, except for two neurotypical control subjects for which we did not receive an SCQ. For one of the neurotypical control subjects without SCQ data, we did collect the Social Responsiveness Scale (SRS) which revealed the subject did not have any autism characteristics. The PW-UPD group shows two distinct clusters, those scoring above 15 (possible ASD) and those within normal limits (Figure 1).

Four DPSC lines representing each of the PWS subgroups (PW-del, PW-UPD -ASD, and PW-UPD + ASD) and 4 neurotypical control subject lines were selected from the repository for RNAseq studies (Supplementary Table 1). Deciduous (“baby”) teeth from neurotypical control subjects were collected locally from the Department of Pediatric Dentistry and Community Oral Health at UTHSC and PWS subjects were collected remotely by parents and guardians of the subjects using our previously published protocols (Goorha and Reiter, 2017). DPSC lines were differentiated into neurons using a published neuronal differentiation protocol (Kiraly et al., 2009; Goorha and Reiter, 2017). DPSC from all groups show positive expression for the stem cell marker, NESTIN, while these same lines were negative for neuronal markers MAP2 and GABA A receptor subunit beta 3 (Figure 2A). After the neuronal differentiation protocol and a 3-week maturation period, the cultures were positive for both MAP2 and GABA A receptor subunit beta 3 (Figure 2A). Dental pulp stem cells show a flat, fibroblast like morphology (Figure 2B), but after a month of differentiation, the cultures display a pyramidal neuron-like morphology (Figure 2C). Both neurotypical and PWS neuronal cultures displayed similar morphologies and cell numbers throughout the differentiation process. We have previously established that a 3-week neuronal maturation period is sufficient to induce significant gene expression changes indicative of terminal neuronal differentiation (Urraca et al., 2015) and reflect underlying disease specific gene expression patterns (Urraca et al., 2018).

RNAseq files from each individual were analyzed by the UTHSC Bioinformatics Core to create lists of significantly (p-value ≤ 0.05 and FDR ≤ 0.05) different transcripts for each PWS subtype versus control and each subtype versus the other genotypes (see methods for details). Using these lists of significantly differentially expressed transcripts for each group versus control, we created Venn diagrams in BioVenn (Hulsen et al., 2008) (Figure 3A). From this analysis we identified a unique molecular signature for each subgroup and a core PWS signature comprised of 3 transcripts within the PWS critical region (SNRPN, SNURF, and MAGEL2) and 9 transcripts coding for genes located outside of 15q PWS/AS critical region (Figure 3B). Supplementary Table 2 lists these PWS specific transcripts, their function, and evidence that ties these genes or their function to PWS phenotypes. Analyzing the expression of genes across the PWS/AS critical region in our subjects shows the expected expression changes from control to PWS subtypes (Supplementary Figure 1). Specifically, expression of genes that are maternally imprinted in PWS (MAGEL2, SNRPN, SNURF, NDN, and MRKN3) are absent in our PWS groups, regardless of genetic subtype. These data confirm that DPSC neurons from subjects accurately recapitulate the genetic landscape of PWS in terms of imprinted gene expression.

In order to find molecular changes specific to ASD in our dataset, we investigated the significantly different transcripts (p-value ≤ 0.05, FDR ≤ 0.05) in the -ASD groups versus the UPD + ASD group (Figure 4A). 380 transcripts met our significance cut off and were interrogated for enrichment analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID), which assigns enrichment scores for Gene Ontology (GO) and other descriptive terms significantly over-represented in our differentially expressed dataset as compared to the entire human genome (Huang da et al., 2009a, b). Using DAVID functional clustering for GO terms, we found significant enrichments in mitochondrial compartments, functions, and processes (p-value ≤ 0.002) (Figure 4B). The top enrichment clusters were mitochondrial membrane transcripts and transcripts involved in mitochondrial maintenance and biogenesis. It is important to note that no other enrichments (score ≥ 3.0) were found in this dataset. From the transcripts identified in the DAVID enrichment analysis, we created a heatmap across all groups for the mitochondrial genes responsible for the enrichment scores using ClustVis (Metsalu and Vilo, 2015) (Figure 4C). This heatmap illustrates that most mitochondrial transcripts identified in our enrichment analysis are dramatically decreased in the PW-UPD + ASD group. Supplementary Table 3 lists these transcripts and their function. These data support our premise that mitochondrial dysfunction may contribute to increased ASD risk in PW-UPD sub-group.

Since enrichment analysis revealed global down regulation of mitochondrial transcripts, we looked for mitochondrial defects at the cellular level in PW-UPD + ASD neurons. Anti-Beta Tubulin was used to outline the area of the neurons (Figures 5A,C). To visualize mitochondria in neurons, we used the mitochondrial surface marker TOMM20 (translocase of outer mitochondrial membrane 20) (Lavie et al., 2017; Nguyen et al., 2018). Anti-TOMM20 revealed a mitochondrial phenotype within the + ASD neurons characterized by perinuclear aggregation of mitochondria and decreased mitochondria detected in the neuronal projections (Figures 5B vs. D). To quantify this mitochondrial specific phenotype, we performed a blinded analysis of neuron images from neurotypical control subjects and PWS subjects. For each group ≥5 individuals were used. Each cell line was coded to ensure no user bias when measuring the mitochondrial area within the neuron. The neuronal cultures were immunolabeled with anti-TOMM20 and anti-Beta Tubulin and then imaged through the cell by confocal microscopy. To determine the percentage of neuronal area containing mitochondria, Imaris software (Oxford Instruments, Abingdon, United Kingdom) was used on the confocal stacks. Neuron images were uploaded to the Imaris program and surfaces for quantification were created using the anti-Beta Tubulin (Figures 6A,C) signal to compute total neuronal area and anti-TOMM20 (Figures 6E,G) to compute the mitochondrial volume within the neuron. Figure 6 shows the Imaris output after using the “surfaces” function to define the area of both components. We were able visually to confine our measurements to single neurons and eliminate any debris from this analysis. For each cell line, ≥15 neurons were imaged and analyzed. We found a significant difference in the average neuronal mitochondrial coverage between the PW-UPD + ASD group and the other groups (Figure 7). These data confirm the hypothesis that mitochondrial dysfunction may underlie the increased ASD incidence of the UPD class.