Bioinformatics analyses show dysregulation of calcium-related genes in Angelman syndrome mouse model

BACKGROUND
Angelman syndrome (AS) is a genetic neurodevelopmental disorder caused by the loss of function of the UBE3A protein in the brain. In a previous study, we showed that activity-dependent calcium dynamics in hippocampal CA1 pyramidal neurons of AS mice is compromised, and its normalization rescues the hippocampal-dependent deficits. Therefore, we expected that the expression profiles of calcium-related genes would be altered in AS mice hippocampi.


METHODS
We analyzed mRNA sequencing data from AS model mice and WT controls in light of the newly published CaGeDB database of calcium-related genes. We validated our results in two independent RNA sequencing datasets from two additional different AS models: first one, a human neuroblastoma cell line where UBE3A expression was knocked down by siRNA, and the second, an iPSC-derived neurons cell line from AS patient and healthy donor control.


FINDINGS
We found signatures of dysregulated calcium-related genes in AS mouse model hippocampus. Additionally, we show that these calcium-related genes function as signatures for AS in other human cellular models of AS, thus strengthening our findings.


INTERPRETATION
Our findings suggest the downstream implications and significance of the compromised calcium signaling in Angelman syndrome. Moreover, since AS share similar features with other autism spectrum disorders, we believe that these findings entail meaningful data and approach for other neurodevelopmental disorders, especially those with known alterations of calcium signaling.


FUNDING
This work was supported by the Angelman Syndrome Foundation and by the Israel Science Foundation, Grant Number 248/20.


Introduction
Angelman syndrome (AS) is a genetic neurodevelopmental disorder, manifested by severe cognitive and motor impairments (Angelman, 2008;Williams et al., 2006;Peters et al., 2004;Williams, 2005;Summers et al., 1995). Its prevalence is estimated to be between 1:10,000-1:20,000 (Williams et al., 2006;Dan, 2009). The cause for AS is the loss of function of UBE3A protein in the brain, usually due to deletion of small portions of the maternal chromosome 15(q11-13) that contains the UBE3A gene (Kishino et al., 1997;Matsuura et al., 1997;Knoll et al., 1989;Gustin et al., 2010). Knockout of this gene in mice recapitulates many features of AS (e.g. motor dysfunction, aberrant behavior and cognitive deficits) making this model an efficient tool for investigating the disease (Jiang Y hui, Armstrong D, Albrecht U, Atkins CM, Noebels JL, Eichele G, et al., 1998;Miura et al., 2002;Van Woerden et al., 2007;Sonzogni et al., 2018). Using this model, we and others have shown that brain regions implicated in AS deficits correlate with aberrant cellular excitability (Sidorov et al., 2018;Rotaru et al., 2018;Wallace et al., 2012;Gu et al., 2019;Judson et al., 2016;Kaphzan et al., 2012;Kaphzan et al., 2011;Kaphzan et al., 2013;Egawa et al., 2012;Wang et al., 2018), and although not always (Rotaru et al., 2018), in most cases the recovery of the excitability aberrations rescued the behavioral deficits (Sidorov et al., 2018;Gu et al., 2019;Kaphzan et al., 2012;Kaphzan et al., 2013). Recently, we also showed that activitydependent calcium dynamics in hippocampal CA1 pyramidal neurons (PNs) of AS mice is compromised, and a manipulation normalizing this decreased calcium dynamics rescues the hippocampal-dependent deficits . Therefore, we hypothesized that this aberrant activity-dependent calcium signaling in the AS hippocampus is correlated with altered mRNA expression of various calcium-related genes. Since calcium signaling entails a pivotal role in neuronal functioning, we further posited that these effects could possibly be of two types. The first expected effect would be a direct effect on genes that their mRNA expression is regulated by calcium signaling. The second predicted effect is that the expression profile of genes that regulate calcium signaling will be altered via controlled homeostatic processes as a feedback to the aberrant calcium signaling. In the herein study, we investigated several publicly available gene expression datasets and carried out systematic bioinformatics analyses of the expression profiles of calcium-related genes taken from CaGeDB database (Hörtenhuber et al., 2017). In the CaGeDB database, calcium-related genes are divided into two categories, calcium-target genes (genes that are regulated by calcium dynamics) and calcium-regulating genes (genes that regulate calcium signaling and might be altered as a homeostatic response to changes in intracellular calcium levels). First, we analyzed the expression profiles of calcium-related genes in an mRNA transcriptome dataset from hippocampi of AS (Angelman Syndrome) and WT (Wild Type) littermates recently generated by us . We found that the expression of calcium-related genes is altered in AS mouse model hippocampus compared to WT controls. Interestingly, in addition to calcium-related transcriptome differences between healthy and AS model mice we found differences in calcium signaling between sexes both in WT and in AS model mice. Using a machine learning approach of Random Forest (RF), a classification algorithm that has been successfully applied in several recent RNA sequencing studies (Wenric and Shemirani, 2018;Ram et al., 2017), we generated a gene expression signature of calcium-target genes and another gene expression signature of calcium regulating genes that can differentiate between AS and WT. To further validate our findings of dysregulated calcium-related genes in AS, we investigated the herein found calcium-dependent signatures of AS in two additional independent publically available mRNA sequencing datasets generated from human cellular models of Angelman syndrome.

Mouse hippocampus RNA sequencing data
We utilized AS model mice RNA sequencing dataset that we recently generated and published . In short, this dataset consists six WT and six AS littermates with three males and three females in each genotype. All mice were 3-4 months old with either maternal deficiency of Ube3a (AS) or wild type (WT) littermates on C57BL/6 background strain. The genotyping for classification of WT (Ube3a p+/ m+ ) and AS (Ube3a p+/m− ) was conducted by PCR procedure using specific primers as previously described (Jiang Y hui, Armstrong D, Albrecht U, Atkins CM, Noebels JL, Eichele G, et al., 1998). All procedures were performed in strict accordance with the University of Haifa regulations and the US National Institutes of Health guidelines (NIH publication number 8023). Following cervical dislocation, brains were quickly extracted on ice-cold cutting solution (in mM): 110 Sucrose, 60 NaCl, 3 KCl, 1.25 NaH2PO4, 28 NaHCO3, 0.5 CaCl2, 7 MgCl2, 5 Glucose. Hippocampi were immediately isolated, and the ventral third of the right hippocampus was snap frozen in liquid nytrogen. Tissue was lysed in lysis buffer using a probe sonicator. Total RNA was isolated using the RNeasy Lipid Tissue Mini Kit, Cat No: 74804 (QIAGEN) according to manufacturer instructions. Isolated RNA concentration and quality was determined by Qubit® quantitation assay using Qubit® 2.0 fluorometer (Invitrogen -life technologies, USA) and Agilent 4200 TapeStation System (Agilent Technologies, USA). Samples were prepared for Illumina poly-A RNA sequencing using NEB Ultra RNA Library Prep Kit for Illumina (NEB#7530) (BioLabs Inc., MA, USA) according to the manufacturer protocol. Libraries were loaded at 8 PM and were sequenced with a 2 × 100 bp PE run on Illumina HiSeq 2500 using a V3 flow cell.

Bioinformatics analyses of mouse hippocampus RNA sequencing data
RNA sequencing data was processed by the Tauber Bioinformatics Research Center parallelized pipeline. The raw reads were cleaned from adaptors and PCR duplicates using Trimmomatic algorithm (Bolger et al., 2014). The cleaned reads were aligned with STAR algorithm (Dobin et al., 2013) to mouse genome assembly GRCM38:mm10 https: //www.ncbi.nlm.nih.gov/assembly/GCF_000001635.20/ with default parameters. The expression levels were quantified by RSEM algorithm (Li and Dewey, 2011) in FPKM units.
Differential expression analysis was performed with DeSeq2 (Love et al., 2014) algorithm on raw read expression matrices. DeSeq2 normalized tables were used for further bioinformatics analyses. Heat maps were generated using the Heatmapper online tool (Babicki et al., 2016) with complete linkage clustering method and Pearson distance measure.

Random Forest Analysis
We utilized RF to identify signature of genes that based on their expression it was possible to classify between AS and WT control. We added an iterative procedure to the standard random Forest R package (Liaw and Wiener, 2002) to iteratively identify genes most frequently chosen for the prediction model. We used the 1000-tree RF iterative procedure with 1000 iterations (Supplementary File 1). Principle component analysis (PCA) was performed with R MASS package (Venables and Ripley, 2002) to further validate the prediction power of the chosen genes.

Publically available RNA sequencing data
In order to validate our found signatures from mouse hippocampi dataset, we utilized two independent publically available datasets. RNA sequencing data generated from human neuroblastoma cell lines (SH-SY5Y) in which UBE3A was knocked down by utilizing siRNA for UBE3A (siUBE3A) and compared to controlled non-targeting siRNA transfected cells (siCTR) (GEO accession: GSE103309) (Lopez et al., 2017). The second dataset was an RNA sequencing dataset generated from patientderived iPSC derived neurons and compared to iPSC derived neurons derived from a healthy donor (Germain et al., 2014) (SRA accession number: SRP044749).
RNA sequencing data was processed by the Tauber Bioinformatics Research Center parallelized pipeline. The raw reads were cleaned from adaptors and PCR duplicates using Trimmomatic algorithm (Bolger et al., 2014). The cleaned reads were aligned with TopHat2 algorithm (Kim et al., 2013) to human genome assembly GRCh38 https://www. ncbi.nlm.nih.gov/assembly/GCF_000001405.26/.
The expression levels were quantified by RSEM algorithm (Li and Dewey, 2011) in FPKM units. The expression tables were transformed to natural logarithmic scale and quantile normalized with threshold five for further analyses.

Results
As aforementioned, activity-dependent calcium signaling is compromised in AS mice hippocampal CA1 pyramidal neurons . Hence, we hypothesized that the expression of genes that either affect or are directly affected by calcium signaling will be altered in AS mice. To test this hypothesis, we utilized a polyA transcriptome sequencing dataset of ventral hippocampi from AS model mice and their WT littermates , and analyzed the expression profiles of calcium signaling related genes that are collected in the CaGeDB database (Hörtenhuber et al., 2017). The CaGeDB database contains comprehensive lists of calcium related genes that are divided into two groups: genes that regulate calcium (calcium-regulating genes) and genes that their expression or function are dependent on calcium signaling (calcium-target genes). Based on these lists, we performed statistical analyses with both the calcium-target and the calciumregulating genes expressed in AS and WT mice hippocampi in order to understand whether their expression is disrupted in AS. Using random forest machine learning approach, we found gene expression signatures of calcium-target and calcium-regulating genes that differentiate between AS and WT control mice. To validate further the significance of these gene signatures, we examined our findings by investigating two publically available independent datasets. The first study utilized an siRNA knockdown of UBE3A in a human neuronal cell line (SH-SY5Y) (Lopez et al., 2017). This study generated an RNA sequencing dataset of UBE3A knocked-down cells (siUBE3A) and control cells transfected with non-targeting siRNA (siCTR). The second study generated a RNA sequencing data from iPSC-derived neurons from an AS patient and a healthy donor, a technical duplicates from each sample (Germain et al., 2014). Both the calcium-target gene signature and the calciumregulating gene signature successfully differentiated between the human AS model cells and the healthy control cells in the two datasets, thus strengthening their validity as signatures of AS.
To confirm the centrality of calcium-related genes in AS, we culled from the entire list of differentially expressed genes, the genes that are listed in the CaGeDB database (Hörtenhuber et al., 2017) for each sex separately. We found that 27 calcium-related genes are differentially expressed in AS female mice compared to WT female mice, and 42 calcium-related genes are differentially expressed in AS male mice compared to WT male mice (Supplementary Table 3A-B). Interestingly, there was little overlap (4 genes) between the genes dysregulated in AS female and those dysregulated in AS male mice (Fig. 2). Sex-dependent transcriptome differences in healthy and diseased brains have recently become an active research ground and several studies report differences in transcriptome signatures in male and female brains Yagi and Galea, 2019;Bundy et al., 2017). While the abovementioned genes were found by comparing sex-matched WT and AS mice, we also performed a complementary comparison between the sexes within each genotype. This comparison found male to female variation in the expression of calcium-related genes both in WT mice (Supplementary Fig. 1A Table 4B). Surprisingly, in the WT mice we found 69 differentially expressed calciumrelated genes in male compared to female hippocampi, while in AS there were only 23 calcium-related genes differentially expressed between males and females. Except for two genes (Hpcal1 and S100a9), none of the sex-dependent differentially expressed genes in the two genotypes (WT and AS) overlapped. The most significant enriched functional clusters of genes differentially expressed in male compared to female WT mice were 'modulation of chemical synaptic transmission', 'regulation of MAPK cascade', 'Nervous system development', and 'Regulation of synaptic GABAergic transmission' (Supplementary Fig. 1B), while in AS mice the most significant enriched functional cluster was 'Regulation of postsynaptic cytosolic calcium ion concentration' (Supplementary Fig. 1D).

Calcium-target genes
Next, to better understand the dysregulation of calcium-related genes in AS, we separated the calcium-target genes and calciumregulating genes as assigned by the CaGeDB database (Hörtenhuber et al., 2017). Calcium-target genes are genes regulated by transient changes in intracellular calcium concentrations. Due to the centrality of calcium-signaling in neuronal functioning we hypothesized that it is less likely that distinct genes will be strongly affected by alterations in calcium dynamics without any additional homeostatic compensations for such severe interruptions. Therefore, we predicted that the aberrant calcium signaling slightly alters the expression of whole sets of calciumtarget genes. To test this, we applied a multi-run random forest (multirun RF) procedure (see Methods) to identify calcium-target classifier genes that, as a group, can classify samples into AS or WT hippocampi. Because of the previously mentioned discrepancies between calcium signaling in female and male hippocampi, we first found calcium-target signatures for female AS mice, followed by calcium-target signatures for male AS mice, and lastly found a signature for combined female and male cohort of AS mice.

Calcium-target signature for female AS model mice
The multi-run RF analysis performed for WT and AS female mice yielded a list of 10 calcium-target classifier genes ( Fig. 3A, Supplementary Table 5). Based on the expression of these genes, we performed a principle component analysis (PCA) which showed a clear separation between the AS and the WT female mouse hippocampi (Fig. 3B). This clear PCA result confirmed that genes chosen by the multi-run RF procedure are indeed good classifiers for AS and WT female mice hippocampi. From the found genes, four were downregulated (Rhot1, Sparc, Cdh2, Nedd4) and 6 were upregulated (Braf, Gpd2, Egfl7, Syt1, Clstn1, Ncald) in AS female mice compared to WT female mice (Fig. 3A).

Calcium-target signature for male AS model mice
The multi-run RF analysis performed on male WT and AS mice yielded a list of 15 calcium-target classifier genes (Fig. 3C,Supplementary Table 6). As mentioned above, these genes did not overlap with the genes found as calcium-target signature for female mice. Based on the expression of these 15 genes, PCA showed a clear separation between the AS and the WT male mice hippocampi (Fig. 3D), again confirming that genes chosen by the multi-run RF procedure are indeed good classifiers of male AS and WT mice hippocampi. From the found genes, 10 were downregulated (Sulf2, Adcy1, Itsn1, Ahcyl1, Rit1, Unc13b, Creld1, Unc13a, Agrn, Sptan1) and 5 were upregulated (Mcfd2, Rbm22, Ppm1f, Slc25a12, Add1) in AS male mice compared to WT male controls A. B.

Calcium-target signature for AS model mice (males + females)
Next, we wished to identify an overall gene expression signature that differentiates between AS and WT mice independent of sex. For this, we combined the female and male samples and performed similar RF procedure choosing genes that predict AS and WT samples. This procedure identified 30 classifier genes (Fig. 3E,Supplementary Table 7). Interestingly, out of these 30 classifier genes, 12 genes were related to 'synapse' according to GO terms (Adcy1, Agrn, Braf, Dst, Grm5, Macf1, Myo6, Nedd4, Pdcd6ip, Rab3gap1, Syt1, Syt11). Next, we performed three distinct PCAs to validate the 30 calcium-target signature genes and how they differentiate between AS and WT samples: mixed females and males groups (AS versus WT) (Fig. 3F), separately on females (Supplementary Fig. 2A) and separately on males ( Supplementary Fig. 2B). Each of the three PCAs show a clear separation between AS and WT mice.
Interestingly, examining the calcium-target signature genes found by each of the three multi-run RF procedures we observed that there was no overlap between the signature genes yielded from the females-only dataset and the signature genes found from the males-only dataset ( Supplementary Fig. 2C).

Calcium-regulating genes
Due to the compromised calcium dynamics in AS mice hippocampi, we predicted that homeostatic processes will take place and genes that regulate calcium dynamics will be affected as well. Therefore, we A. Venn diagram of downregulated calcium-related (CaGeDB) genes in AS female and AS male mice compared to their WT sex-match controls. Below the diagram, a table indicating the names of genes downregulated in female AS compared to female WT and in male AS compared to male WT. One gene (Calr) was found to be downregulated in both female and male AS compared to their sex-matched controls. B. Venn diagram of dysregulated calcium-related (CaGeDB) genes in AS female and AS male mice compared to their WT sex-match controls. Below the diagram, a table indicating the names of genes upregulated in female AS compared to female WT and in male AS compared to male WT. Three genes (Capn5, Myh7, and Otof)) were found to be upregulated in both female and male AS compared to their sex-matched controls.
Heat map of calcium-target signature genes (10) idenƟfied by random forest procedure in female AS and female WT hippocampi (caption on next page) analyzed the expression profiles of genes that are upstream to calcium signaling and regulate the intracellular Ca 2+ concentrations. For that, we extracted the list of calcium-regulating genes from the CaGeDB database (Hörtenhuber et al., 2017) and utilized the multi-run RF approach on the expression of these genes in the mouse hippocampi dataset. Similar to the above calcium-target analyses, we separately identified the calcium-regulating female signature, then calciumregulating male signature, and the combined males and females signature of calcium-regulating genes that differentiate AS from WT samples.

Calcium-regulating signature for female AS model mice
The multi-run RF analysis performed for WT and AS female mice yielded a list of 10 calcium-regulating classifier genes (Fig. 4A,Supplementary Table 8). Based on the expression of these genes, we performed a PCA, which showed a clear separation between the AS and the WT female mouse hippocampi (Fig. 4B). This clear PCA result confirmed that calcium-regulating genes chosen by the multi-run RF procedure are indeed good classifiers for female AS and WT mouse hippocampi. From the found genes, five were downregulated (Tpt1, Opa1, Gsto1, Hsp90b1, Herpud1) and five were upregulated (Cacna1g, Adcy3, Atp2b2, Ppp3cb, Grin2b) in AS female mice compared to WT female mice (Fig. 4A).

Calcium-regulating signature for male AS model mice
The multi-run RF analysis performed for male WT and AS mice yielded a list of 20 calcium-regulating classifier genes (Fig. 4C,Supplementary Table 9). Interestingly, three genes overlapped with calciumregulating signature genes for female AS mice. Herpud1 gene was downregulated both in AS female and male mice. However, Atp2b2 and Grin2b were upregulated in AS female mice, while downregulated in male AS mice compared to their sex-matched WT controls. PCA based on the expression of the 20 calcium-regulating signature genes showed a clear separation between the AS and the WT male mice hippocampi (Fig. 4D). This further confirms that genes chosen by the multi-run RF procedure are indeed good classifiers for AS genotype. From these 20 found genes, 12 were downregulated and 8 were upregulated in AS male mice compared to WT male controls (Fig. 4C).

Calcium-regulating signature for AS model mice (males + females)
Next, similar as with the calcium-target genes, we wished to identify an overall calcium-regulating signature for the AS genotype independent of sex. For this, we combined the female and male samples and performed similar RF procedure choosing calcium-regulating genes that separate AS and WT samples. This procedure yielded a signature of 30 calcium-regulating classifier genes (Fig. 4E,Supplementary Table 10). PCA based on the expression of these 30 identified calcium-regulating signature genes showed a good separation of AS samples from the WT samples (Fig. 4F). These genes also showed a clear separation for females only (Supplementary Fig. 3A) and males only samples (Supplementary Fig. 3B). Interestingly, from the 30 calcium-regulating signature genes for AS mice, 13 were associated with 'synapse' according to GO terms (Adcy3, Arrb2, Bnip3, Cacna1g, Calm1, Dlg4, Dmd, Fyn, Nptn, Nrxn1, Ppp3cb, Ppp3r1, and Stoml2).
Interestingly, examining the signature genes found by each of the three multi-run RF procedures we observed that only one gene (Her-pud1) was overlapping between the calcium-regulating signature genes yielded independent of sex, females-only dataset, and the signature genes found from the males-only dataset ( Supplementary Fig. 3C).

Validating calcium-target and calcium-regulating signatures: randomly chosen genes
We further examined whether the chosen calcium-target and calcium-regulating signature gene lists (Supplementary Table 7 and  Supplementary Table 10) perform better as signatures of AS than randomly chosen calcium-target and calcium-regulating genes. Utilizing the AS and WT mouse hippocampi dataset we performed PCAs based on the expression of 10 random lists of 30 calcium-related genes. These 10 randomly selected lists of calcium-related genes, 5 lists of 30 calciumtarget genes and 5 lists of 30 calcium-regulating genes, were not able to differentiate between AS and WT ( Supplementary Fig. 4), thus validating the significance of the RF-found calcium-target and calciumregulating gene signatures in mouse hippocampi.

Validating calcium-target and calcium-regulating signatures in human cell models of AS
To validate the significance of the found calcium-target and calciumregulating signatures found independent of sex of the samples, we utilized two independent publically available RNA sequencing datasets. The first dataset that was used for validation is an mRNA sequencing dataset of human neuroblastoma cell lines in which either UBE3A was knocked-down using an siRNA (siUBE3A cells) or a non-target siRNA was used to generate control cells (siCTR) (Lopez et al., 2017). The second dataset that we used for validating the significance of the found signatures was an mRNA sequencing dataset generated from iPSCderived neurons from an AS patient and a matched healthy donor, with two technical replicates per each individual (Germain et al., 2014).

Validating calcium-target signature in human AS models
To examine whether the calcium-target gene signature identified in our AS mouse hippocampi dataset when using mixed samples (both sexes per genotype), can also serve as a signature in other AS models we performed PCA using the expression of these genes in the human cellular model of AS (siRNA in neuroblastoma). Alignment of sequencing reads to the human genome revealed that from the 30 calcium-target signature genes, 28 were expressed in this dataset (Fig. 5A, Supplementary   Fig. 3. Calcium-target gene signatures identified by random forest procedure. A. Heat map of calcium-target signature genes (10) identified by random forest procedure in female AS and female WT hippocampi. Red color indicates downregulated genes and green color indicates upregulated genes. B. PCA of female AS and female WT mouse hippocampi based on expression of calcium-target signature genes (10) identified by random forest procedure in female AS and female WT mouse hippocampi. X and Y-axis correspond to the first two principle components. Red dots indicate AS female mice, green dots indicate WT female mice. C. Heat map of calcium-target signature genes (15) identified by random forest procedure in male AS and male WT hippocampi. Red color indicates downregulated genes and green color indicates upregulated genes. D. PCA of male AS and male WT hippocampi based on expression of calcium-target signature genes (15) identified by random forest procedure in male AS and male WT mouse hippocampi. X and Y-axis correspond to the first two principle components. Red dots indicate AS male mice, green dots indicate WT male mice. E. Heat map of calcium-target signature genes (30) identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. F. PCA of female and male AS and WT hippocampi based on the expression of calcium-target signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle components. Red dots indicate AS mice and green dots indicate WT mice.

WT Male AS Male
Heat map of calcium-regulaƟng signature genes (10) idenƟfied by random forest procedure in female AS and female WT hippocampi A. C.

A. B.
(caption on next page) Table 11). PCA showed that these 28 genes clearly separate the UBE3A knockdown cells from control cells (Fig. 5B) validating the use of these calcium-target genes as a signature of AS. Next, we examined whether the same calcium-target gene signature genes found in our mice dataset, can also separate between iPSC-derived neurons from an AS patient to those of a healthy donor. For this, we aligned the transcriptome reads generated from iPSC-derived neurons and found that all 30 calcium-target signature genes are expressed in this dataset (Fig. 5C, Supplementary Table 12). PCA showed that these genes clearly separate the AS patient samples from the healthy donor ones also in this dataset (Fig. 5D).

Validating calcium-regulating signature in human AS models
Similar to the above validation of calcium-target signature genes we next validated the calcium-regulating signature genes found independent of mice sex. Out of the 30 calcium-regulating signature genes that were identified in the mice dataset, 28 were found in the human neuroblastoma cell line (Fig. 6A, Supplementary Table 13). PCA of these 28 genes showed a very good separation between the UBE3A knockdown cells and the control cells (Fig. 6B).
Likewise, out of the 30 calcium-regulating signature genes that were identified in our mice data, 28 genes were found in the iPSC-derived neurons dataset (Fig. 6C, Supplementary Table 14). Interestingly, these were the same 28 genes that were found in the neuroblastoma cell line. PCA of these genes showed a clear separation between the iPSCderived neurons of AS patient to those of the healthy control donor (Fig. 6D). Again, these analyses validated the use of these calciumregulating genes as a signature of AS.
Interestingly, in spite of these genes being a valid and reliable signature for AS, the direction of dysregulation in the different models was not always the same. For example, some of the signature genes that were found to be upregulated in the samples of AS mice were downregulated in the neuroblastoma-siRNA model. Nonetheless, in both models these genes indicated well the the AS condition separating it from the control samples. The same is true between the human models; some signature genes that were upregulated in the AS iPSC-derived neurons were downregulated in the siUBE3A neuroblastoma cells that model AS.

Protein-protein interaction of found signature genes
To better delineate the relations between all of the 60 calcium-target and calcium-regulating signature genes we utilized the Search Tool for the Retrieval of Interacting Genes (STRING) online tool (Franceschini et al., 2013) to build a protein-protein interacting (PPI) network. The network consisted of 60 nodes and 87 edges (Fig. 7). Seven genes with ≥5 edges were defined as hub genes (Calm1, Dlg4, Hsp90b1, Hspa5, Braf, Syt1, and Ppp3cb) (Supplementary Table 15). The genes with the highest number of edges were Calm1 with 13 edges and Dlg4 with 10 edges. The functions of these two genes are closely related to glutamatergic synaptic plasticity, which is aberrant in AS hippocampus (Jiang Y hui, Armstrong D, Albrecht U, Atkins CM, Noebels JL, Eichele G, et al., 1998;Weeber et al., 2003).

Discussion
In a previous study, we showed that calcium dynamics is compromised in the hippocampi of AS model mice . This reduced activity-dependent calcium signaling was shown to result from the enhanced expression and activity of the α1-Na + /K + -ATPase (α1-NaKA). Normally, the α1-NaKA activity in the hippocampus is marginal, while the α3-NaKA isoform is dominant (Chakraborty et al., 2017). However, in AS hippocampus the α3 expression is normal but the α1 isoform is significantly enhanced (Kaphzan et al., 2013). We showed that inhibiting the increased activity of α1-NaKA in AS mice restored calcium dynamics and rescued the hippocampal dependent deficits, probably via slowing the sodium gradient used for evacuating the activity-dependent calcium influx . Hence, we predicted that these differences in calcium dynamics are reflected in altered expression of calcium-dependent genes. To test this hypothesis we analyzed a mouse hippocampal mRNA sequencing dataset of AS mice and their WT littermates  in light of a publicly available database of calcium-related genes (CaGeDB) (Hörtenhuber et al., 2017).
Initially we took an unsupervised approach and found genes that are differentially expressed in AS male mice and AS female mice compared to sex-matched WT controls. Enrichment cluster analysis of these genes showed two highly significant clusters that are directly related to calcium signaling (Fig. 1). Next, we focused on calcium-related genes from a curated database that includes all calcium-signaling-related genes, CaGeDB (Hörtenhuber et al., 2017). Again, we examined the expression profiles of genes in male and female AS and WT mice. Many genes were differentially expressed both in female AS mice compared to female WT mice and in male AS mice compared to male WT mice, however, there was little overlap between the differentially expressed genes in female AS mice and those in the male AS mice compared to their sex-matched controls (Fig. 2). Because of this discrepancy between females and males, we further examined the female to male variation within each genotype. Interestingly, we found that not only the female to male variation in WT was different from the variation in AS, but also that the number of genes that vary between the sexes in WT were 3 times bigger than the number of genes that vary in the AS. These transcriptomic differences in female to male variations coincide with our previous findings regarding behavioral phenotypes in AS mice . These findings showed that some behaviors, such as right-sided exploration, social odor preference and odor sensitivity, have a substantial male to female variation in WT mice, while AS mice do not exhibit male to female variation. Interestingly, the opposite was not observed; there were no behavioral male to female variations in AS that Fig. 4. Calcium-regulating gene signatures identified by random forest procedure A. Heat map of calcium-regulating signature genes (10) identified by random forest procedure in female AS and female WT hippocampi. Red color indicates downregulated genes and green color indicates upregulated genes. B. PCA of female AS and female WT mouse hippocampi based on expression of calcium-regulating signature genes (10) identified by random forest procedure in female AS and female WT mouse hippocampi. X and Y axis correspond to the first two principle components. Red dots indicate AS female mice, green dots indicate WT female mice. C. Heat map of calcium-regulating signature genes (20) identified by random forest procedure in male AS and male WT hippocampi. Red color indicates downregulated genes and green color indicates upregulated genes. D. PCA of male AS and male WT hippocampi based on expression of calcium-regulating signature genes (20) identified by random forest procedure in male AS and male WT mouse hippocampi. X and Y-axis correspond to the first two principle components. Red dots indicate AS male mice, green dots indicate WT male mice. E. Heat map of calcium-regulating signature genes (30) identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. F. PCA of both female and male AS and WT hippocampi based on expression of calcium-regulating signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle components. Red dots indicate AS mice and green dots indicate WT mice. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) A. B.

C. D.
Heat map of calcium-target signature genes (30) expressed in human neuroblastoma cell line treated with siUBE3A (AS model) or non-specific siRNA (control) PCA of siUBE3A (AS model) and non-specific siRNA (control) based on expression of calcium-target signature genes (30) idenƟfied by random forest procedure in both female and male AS and WT mouse hippocampi PCA oĬuman iPSC-derived ne urons from AS paƟent and from healthy control based on expression of calcium-target signature genes (30) idenƟfied by random forest procedure in both female and male AS and WT mouse hippocampi Heat map of calcium-target signature genes (30) expressed in human iPSC-derived neurons from AS paƟent and from healthy control . The calcium-target signature genes (30) were identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. B. PCA of siUBE3A (AS model) and non-specific siRNA (control) based on expression of calcium-target signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle. C. Heat map of calcium-target signature genes (30) expressed in human iPSC-derived neurons from AS patient and from healthy control. The calcium-target signature genes (30) were identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. D. PCA of human iPSC-derived neurons from AS patient and from healthy control based on expression of calcium-target signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) were absent in WT mice. Taken together, these findings show that sex differences play a role beyond genotype differences with regard to calcium-related genes.
In the CaGeDB database, the calcium-related genes are divided into two major types, genes that are downstream targets of calcium signaling and are affected by transient changes in intracellular calcium concentrations (calcium-target genes) and genes that are upstream to calcium signaling and regulate the transient changes in intracellular calcium concentrations (calcium-regulating genes).
The approach of multi-run RF successfully identified sets of genes A. B.

C. D.
Heat map of calcium-regulaƟng signature genes (30) expressed in human neuroblastoma cell line treated with siUBE3A (AS model) or non-specific siRNA (control) PCA of siUBE3A (AS model) and non-specific siRNA (control) based on expression of calcium-regulaƟng signature genes (30) idenƟfied by random forest procedure in both female and male AS and WT mouse hippocampi  (30) were identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. B. PCA of siUBE3A (AS model) and non-specific siRNA (control) based on expression of calcium-regulating signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle. C. Heat map of calcium-regulating signature genes (30) expressed in human iPSC-derived neurons from AS patient and from healthy control. The calcium-regulating signature genes (30) were identified by random forest procedure in both female and male AS and WT mice. Red color indicates downregulated genes and green color indicates upregulated genes. D. PCA of human iPSC-derived neurons from AS patient and from healthy control based on expression of calcium-regulating signature genes (30) identified by random forest procedure in both female and male AS and WT mouse hippocampi. X and Y-axis correspond to the first two principle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) that, as a group, differentiate between the two genotypes. These results suggest that sometimes the differential expression of a single gene does not entail statistical significance on its own, but as an aggregate of genes, these small incremental changes in expression make a difference. Using RF method, we extracted from the calcium-target genes a list of 30 classifiers that based on their expression profiles it is possible to differentiate clearly between the WT and AS genotypes regardless of sex ( Fig. 3E-H, Supplementary Table 7). Moreover, from the calciumregulating genes we extracted another list of 30 classifiers that can also differentiate clearly between the WT and AS genotypes regardless of sex ( Fig. 4E-H, Supplementary Table 10).
To show that these genes lists are valid as signatures of AS, we also examined whether any random lists of 30 calcium-related genes can function as a signature of AS genotype. We generated 10 random lists of calcium-related genes, 5 lists of 30 calcium-target genes and 5 lists of 30 calcium-regulating genes. Using the AS and WT mouse hippocampi dataset we performed PCAs based on the expression of these randomly chosen lists of calcium-target and calcium-regulating genes. These randomly selected calcium-target and calcium-regulating genes lists were not able to differentiate between AS and WT ( Supplementary  Fig. 4), thus validating the significance of the RF-found calcium-target and calcium-regulating gene signatures.
To further validate the chosen calcium-target and calcium-regulating genes as signatures for AS, we utilized two publically available datasets of human cellular models of AS (Lopez et al., 2017;Germain et al., 2014). We found that both of the signature genes lists (target and regulating), found in the AS mouse model, function well as signatures of AS in these human AS cellular models (Figs. 5-6). The fact that these genes function as a signature in other AS models further emphasize their involvement in AS. The full meaning for the involvement of these specific genes in AS is not yet known. Nonetheless, these genes might be utilized for further understanding of the entire pathological alterations induced by UBE3A deletion. From the herein results, and from additional studies we performed Panov et al., 2020), it appears that UBE3A regulates many processes, some of which are calcium dependent. For example, neuronal spine formation (Cornelia Koeberle et al., 2017), a process highly regulated by calcium , is known to be affected in AS, as observed by the reduced spine density in neurons of AS mice (Dindot et al., 2008). Furthermore, there are many endeavors by us and other groups aimed at rescuing the AS deficits (Van Woerden et al., 2007;Kaphzan et al., 2012;Rayi et al., 2019;Meng et al., 2015;Sun et al., 2015;Huang et al., 2012;Meng et al., 2013). However, it is difficult to assess the efficacy of these interventions before their application in humans. We propose that the abovementioned gene signatures may be seen as endophenotypes that can serve as assessment tools or markers for the efficiency of any therapeutic strategy aiming to rescue the AS phenotype, whether it is in a mouse model or in any other cellular model of AS.
Protein-protein interacting network revealed several hub genes, with the most prominent being Calm1 and Dlg4 (Fig. 7). Interestingly, several The network nodes indicate the calcium signature genes identified by random forest procedure. Differently colored lines represent different types of evidence used in predicting the gene associations. A green line indicates a neighborhood evidence; a blue lineliterature co-occurrence evidence; a purple line -experimental evidence; a light blue line -database evidence. The network consisted of 60 nodes and 87 edges. Seven genes with ≥5 edges were defined as hub genes. Red stars indicate hub genes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) studies showed that Calm1 directly interacts with Ube3a (Jung et al., 2005;Martinez-Noel et al., 2012;Kapoor et al., 2016). Calm1 encodes for the calmodulin protein, which is an intermediate calcium-binding messenger. Activation of neurons induces calcium ions influx. These calcium ions bind to calmodulin generating an active calcium/calmodulin protein, which further binds proteins and activates them (Wayman et al., 2008). In neurons specifically, calcium/calmodulin activates one of the main proteins required for learning and memory, the calcium/ calmodulin-dependent protein kinase II (CaMKII) (Wayman et al., 2008). CaMKII binds to the postsynaptic density (PSD), where it is active and affects synaptic plasticity. The major PSD scaffold protein is PSD95, which is encoded by Dlg4 gene (Gardoni et al., 2006;Li et al., 2016;Li et al., 2017), the second prominent hub protein. Interestingly, a study of AS mice hippocampal pathophysiology tied these two proteins together by showing a decreased activity of CaMKII together with a significant reduction in the association of CaMKII to PSD (Weeber et al., 2003).
To conclude, the herein study was performed in order to further elucidate the molecular changes associated with our previous finding of a compromised calcium signaling in AS mice hippocampi. For that, we examined whether the expression of calcium-related genes is altered in AS. By utilizing unsupervised and supervised approaches, we showed that the expression profiles of distinct calcium-related genes can be used as signatures of AS. Furthermore, we showed that these distinct genes are also reliable classifiers that can clearly distinguish between AS and WT in two independent publicly available datasets of human cellular models for Angelman syndrome. This validates these distinct genes as signatures of AS.
In addition to the abovementioned findings, our approach of emphasizing the centrality of calcium signaling in AS, might be relevant for other neurodevelopmental disorders beyond AS, for it is well established that other autism spectrum disorders are strongly related to dysregulated calcium signaling (Nguyen et al., 2018;Schmunk et al., 2017). We suggest that manipulating calcium signaling may serve as a strategy in the armamentarium of therapeutic interventions for neurodevelopmental disorders. Moreover, transcriptomic signatures of calcium-related genes might serve as markers for the success of therapeutic interventions when tested in neurodevelopmental disease models.

Availability of data
No new data was generated for this study. The analysis was performed on data from the following SRA datasets: Mouse hippocampi data: https://www.ncbi.nlm.nih.gov/bioproject /PRJNA484226 Human neuroblastoma siRNA data: https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE103309 Human iPSC data: https://trace.ddbj.nig.ac.jp/DRASearch/study? acc=SRP044749 Author contributions HK initiated the study; JP and HK were in charge of the concepts and ideas of the research, study design, analysis of data, interpretation of analysis and writing the manuscript. HK was in charge of getting the resources for the work.

Declaration of Competing Interest
All of the authors declare having no conflicts of interests, and are not related to any financial organization or patents.