RNA sequencing and proteomics approaches reveal novel deficits in the cortex of Mecp2-deficient mice, a model for Rett syndrome

Animals

All experimental protocols were followed according to NIH guidelines and approval from the Animal Care and Use Committee of the University of Alabama at Birmingham. Wild-type (WT) males were bred with heterozygous Mecp2 ^tm1.1Jae/+ (Jaenisch) female mice [40]. The founding heterozygous Mecp2 ^tm1.1Jae/+ female mice and WT male mice (C57BL/6J) were obtained from the Mutant Mouse Resource and Research Centers (MMRRC). Colonies were occasionally refreshed with WT male mice from Charles River Laboratories (C57BL/6NCrl) for breeding with Mecp2 ^tm1.1Jae/+ female mice. Genotypes of the resulting offspring were validated by PCR of DNA collected from tail clips. Mutant male mice (Mecp2 ^Jae/y ) along with their WT littermates were collected after postnatal day 60 (P60⁺), which represents the time point when Mecp2-deficient males are considered to be “symptomatic.” We define symptomatic as exhibiting hind limb clasping upon suspension from tail and decreased mobility [40], both of which were observed in all Mecp2-deficient males at the time of collection. For both RNA sequencing and proteomics experiments (methods outlined below), an n of 4 biological replicates per genotype were used, with WT littermate animals used as controls. The number of biological replicates was chosen based on a previous study examining the relationship between the number of biological replicates, RNA-sequencing depth, and statistical power [42].

RNA and protein isolation

Animals were anesthetized under CO₂ followed by rapid decapitation. Whole brains were dissected in ice-cold phosphate-buffered saline and then separated into the cerebellum, brain stem, midbrain, and hippocampus. One hemisphere of the whole cortex was dedicated for RNA and the other for protein isolation. RNA was isolated by placing the cortical hemisphere dedicated for RNA in 1 mL of RNAlater™ solution (Invitrogen/Thermo Fisher Scientific) and allowed to react at 4 °C for at least 1 week. RNA was collected using the Invitrogen™ Ambion™ PureLink™ RNA Mini Kit (Fisher Scientific) according to the manufacturer’s instructions, with the following modifications: (1) tissue was homogenized in 1 mL lysis buffer with β-mercaptoethanol (Sigma-Aldrich) in a dounce homogenizer 5 times and rested on ice for 5 min, followed by another 5 rounds of homogenization; (2) approximately 15 mg of the homogenate was removed and brought up to a final volume of 600 μL in lysis buffer with β-mercaptoethanol (Sigma-Aldrich), then re-homogenized as described above. RNA was eluted in 30 μL of autoclaved and filtered Mill-Q® water. Proteins were isolated from the remaining cortical hemisphere in ice-cold lysis buffer (100 mM Tris base, pH 7.5 at room temperature, 1% (w/v) SDS) with protease inhibitor cocktail and phosphatase inhibitor cocktail 3 (Sigma-Aldrich, product numbers P8340 and P0044, respectively) such that the final concentration was 40 mg/mL. Lysates were sonicated using the Model 120 Sonic Dismembrator (Fisher Scientific) for 7 s at 70% amplitude, pulse 20 s, and rest 50 s for 2 cycles. Protein lysates were centrifuged to pellet debris. The supernatant was removed into another tube and quantified using a bicinchoninic acid (BCA) assay kit (Pierce/Thermo Scientific).

RNA sequencing

RNA samples were submitted to the Genomics Core Laboratory in the Heflin Center of Genomic Sciences at the University of Alabama at Birmingham for sample preparation and sequencing. The samples were first DNase-treated and assessed for total RNA quality using the Agilent 2100 Bioanalyzer, followed by 2 rounds of polyadenylate positive (poly A+) selection and conversion to cDNA. RNA sequencing was performed on the Illumina HiSeq 2500 using the latest versions of sequencing reagents and flow cells, providing up to 300 GB of sequence information per flow cell. TruSeq library generation kits were used according to the manufacturer’s instructions (Illumina). Library construction consisted of random fragmentation of the poly A+ mRNA, followed by cDNA production using random primers. The ends of the cDNA were repaired, A-tailed, and adaptors ligated for indexing (up to 12 different barcodes per lane) during the sequencing runs. The cDNA libraries were quantitated using qPCR in a Roche LightCycler 480 with the Kapa Biosystems kit for library quantitation (Kapa Biosystems) prior to cluster generation. Clusters were generated to yield approximately 725K–825K clusters/mm². Cluster density and quality were determined during the run after the first base addition parameters were assessed. Paired-end 2 × 50 bp sequencing runs were performed to align the cDNA sequences to the reference genome mouse mm10. Approximately 15 million paired 50 bp reads were obtained per sample.

RNA-Seq bioinformatics analysis

All RNA-Seq fastq reads (GEO Series accession number GSE96684) were processed in the Galaxy platform [43]. First, raw RNA-Seq reads were concatenated using the “Concatenate datasets tail-to-head” tool. The concatenated raw fastq reads were then trimmed using Trim Galore! (Galaxy version 0.4.2; [44]) to remove adapter sequences and low-quality base pairs with the following parameters: selected “paired-end” library; trimming reads—automatic detection; trims 1 bp off every read from its 3′ end—yes; no removal of N bp from 3′ end of reads 1 and 2 (respectively); and all other settings left at default. Trimmed fastq reads were then run through FastQC (Galaxy version 0.65; [45]) for additional quality control measures. Following quality control, the reads were aligned to the mouse mm10 reference genome using TopHat (Galaxy version 2.1.0; [46, 47]) using the following parameters: mean inner distance between mate pairs—175; standard deviation for distance between mate pairs—20; report discordant pair alignments—yes; and all other settings left at default. Aligned reads were assembled using Cufflinks (Galaxy version 0.0.7; [46, 48]) using the following parameters: max intron length—300,000; min isoform fraction—0.1; pre-mRNA fraction—0.15; perform quartile normalization—yes; use reference annotation—use reference annotation, with the iGenomes UCSC mm10 genome used as the reference annotation [49]; perform bias correction—yes, reference sequence data—locally cached, using reference genome—mouse mm10; use multi-read correct—yes; and job resource parameters—left at default values. Following transcript assembly and estimated fragments per kilobase of transcript per million fragments mapped (FPKM) abundances, each sample’s assembled transcript was merged using the Cuffmerge (Galaxy version 2.2.1.0; [46, 48]) tool with the following parameters: use reference annotation—yes, reference annotation—iGenomes UCSC mm10 genome [49]; use sequence data—no; minimum isoform fraction—0.05; and job resources parameters—left at default. Finally, Cuffdiff (Galaxy version 2.2.1.3; [46, 48]) was used to calculate statistical changes in gene expression using the following parameters: transcripts—output file from Cuffmerge step; omit tabular data sets—no; generate SQLite—yes; input data type—SAM/BAM; condition 1—WT TopHat accepted hits files; condition 2—Mecp2 ^Jae/y TopHat accepted hits files; library normalization method—quartile; dispersion estimation method—per-condition; false discovery rate—0.05; minimum alignment count—100; use multi-read correct—yes; perform bias correction—yes; reference sequence data—locally cached, reference genome mouse mm10; include read group data sets—yes; include count based output files—yes; apply length correction—cufflinks effective length correction; and all other remaining parameters were left at default settings. Since poly A+ selection was utilized to generate the cDNA libraries used for RNA-Seq, any significant and differentially expressed genes that mapped to a putative non-coding gene were removed from analysis.

Proteomics

Western blotting

Lysates were boiled in 2× sample loading buffer (100 mM Tris base, pH 6.8, 4% SDS in Laemmli-sodium dodecyl sulfate, 600 mM β-mercaptoethanol, 200 mM dithiothreitol (DTT), and 20% glycerol) in a 1:1 ratio for 15 min at 60 °C. A final concentration of 10 μg per sample was loaded into a 4–20% gradient pre-cast mini-PROTEAN® TGX™ gel (Bio-Rad) and ran at 200 V in 1× running buffer (24.76 mM Tris base, 190 mM glycine, 0.1% SDS).

For blots tested for RNPEP, QDPR, and CIRBP protein expression, gels were transferred to a nitrocellulose membrane using the Trans-blot turbo system (Bio-Rad), mixed molecular weight protocol (2.5 A, 25 V for 7 min), followed by blocking with LI-COR® blocking buffer at a 1:1 ratio with TBS. The blots were incubated with primary polyclonal chicken anti-GAPDH (EMD Millipore) at 1:2500 for 15 min at room temperature, while the remaining primary antibodies were incubated overnight at 4 °C: polyclonal rabbit anti-RNPEP (Proteintech) at 1:1000, polyclonal rabbit anti-QDPR (Proteintech) at 1:1000, and polyclonal rabbit anti-CIRBP (Proteintech). All secondary antibodies were either goat anti-rabbit or goat anti-chicken (LI-COR®), all incubated at 1:10,000 for 1 h at room temperature. Imaging was performed on a LI-COR® Odyssey machine with 1 or 1.5 intensity on both the 680 and 800 channels. All rabbit antibodies were imaged in the 680 channels and chicken in 800.

For blots tested for MeCP2, NEFM, mGluR3, and HSPH1 protein expression, the gels were transferred to an Immobilon-P PVDF membrane (EMD Millipore) for 1 h at 100 V in 1× transfer buffer (250.7 mM Tris base, 190 mM glycine, 20% methanol) with an ice pack at room temperature, followed by blocking with 10% non-fat dry milk/0.1% Tween-20/TBS (10% milk-TBST) for 1 h at room temperature. The blots were initially incubated with primary antibody chicken anti-GAPDH at 1:2000 (EMD Millipore) for 1 h at room temperature, washed 3 times in 10% milk-TBST, and incubated for 1 h at room temperature with secondary antibody goat anti-chicken at 1:2000 (Santa Cruz Biotechnology). Blots were then washed 3 times with 10% milk-TBST and imaged with Classico chemiluminescent reagent (EMD Millipore) using an autoradiography film developing system (Denville Scientific). The blots were then individually probed with the following antibodies and conditions: mouse anti-MeCP2 at 1:2000 (Sigma) for 2 h at room temperature, rabbit anti-NEFM at 1:1000 (Proteintech) for 2 h at room temperature, rabbit anti-mGluR3 at 1:500 (Alomone labs) overnight at 4 °C, and rabbit anti-HSPH1 at 1:150 (Novus Biologicals) overnight at 4 °C. The respective secondary antibodies were incubated with blots at 1:2000 for 1 h at room temperature: goat anti-mouse (Santa Cruz) and goat anti-rabbit (Santa Cruz). Protein quantification for all blots was performed using Image Studio Lite (LI-COR®) and Origin2015 (OriginLab). Relative protein amounts per antibody tested were normalized to GAPDH expression in the same lane.

Pathway analysis

Ingenuity® Pathway Analysis (IPA®; Qiagen) was used to identify significant biological pathways in both RNA-Seq and proteomics data sets. A list of detected genes and detected proteins (including post-translational modifications) was used as the data input for both individual and comparison pathway analyses, using a q < 0.05 cutoff for the gene pathway and p < 0.1 cutoff for the protein pathway analyses [54] such that only significant genes/proteins were considered for significant pathways. The “User dataset” option was chosen to use each individual detected gene/protein data set as the “reference set” for which to generate significant pathways. Pathways from the “diseases and biological functions” category were used for comparison analyses. Fisher’s t-test of p < 0.05 (or −Log₁₀ p-value > 1.3) was used to determine statistical significance of a pathway.

Global gene expression in symptomatic Mecp2 ^Jae/y whole cortex

We began our studies by performing RNA-Seq in whole cortical tissue obtained from the Mecp2 ^Jae/y murine model of RTT. Previous studies have analyzed the transcriptome of RTT mice; however, analysis performed in the same affected brain region in the same animal model allowed us to make general comparisons across RNA-Seq and proteomics platforms. While each data set offers unique insight into the disease, the comparison of the two allows us to make general assumptions that (1) differentially expressed “hits” common to both platforms possibly share a transcriptional mechanism of dysregulation and (2) proteins differentially expressed without a similar change in gene expression may indicate differences in either protein stability or posttranscriptional regulation. We chose to analyze whole cortex (WCX) in a Mecp2 ^Jae/y murine model of RTT because MeCP2 protein is highly expressed in WCX [10, 11] and has pathological characteristics of RTT [26]. RNA sequencing (RNA-Seq) was used to quantify gene expression of 15 million paired 50 bp reads in WT and symptomatic Mecp2 ^Jae/y mice (as described in the “Methods” section). We identified 391 significant, differentially expressed (DE) genes (Fig. 1a, Additional file 1). Of these 391 genes, 132 genes had increased expression and 259 had decreased expression. To assess the consistency of our data with previously published microarray and RNA-Seq studies of RTT, we compared our findings to transcriptome-based data on different species and brain regions [4, 19–26, 33, 34, 55, 56]. In our data, 35 genes were identified as known RTT differentially expressed genes and are represented in Fig. 1b and Table 1. Of the 35 genes, 18 were identified in supplemental material from [4] and [34] (Table 1). Additionally, within the group of 35 identified RTT hits, we found that 66% of the DE genes had decreased expression, while 34% had increased expression (Fig. 1b, Table 1). These observations are consistent with previous studies with the exception of one gene Cacnb3. In our data set, Cacnb3 is significantly decreased; however, Tudor et al. found that this same gene was upregulated [20].

Gene targets: this study	UniProt: function keywords	Fold change	RTT gene hits: direction of expression (with references)
Aacs*	• Molecular function: ligase • Biological process: fatty acid metabolism; lipid metabolism • Ligand: ATP-binding; nucleotide-binding	− 0.79	Decrease4
Aldh1a1*	• Molecular function: oxidoreductase • Ligand: NAD	− 0.67	Decrease4
Arhgdig	• Molecular function: GTPase activation	− 0.69	Decrease 2
Auts2	N/A	1.40	Increase7
C1qa	• Biological process: complement pathways; immunity; innate immunity	− 0.76	Decrease7
Cacnb3	• Molecular function: calcium channel; ion channel; voltage-gated channel • Biological process: calcium transport; ion transport; transport • Ligand: calcium	− 0.79	Increase 2
Calb1	• Ligand: calcium; metal-binding; vitamin D	− 0.79	Decrease 5
Calr*	• Molecular function: chaperone • Ligand: calcium; lectin; metal-binding	− 0.80	Decrease8
Dgkg*	• Molecular function: kinase; transferase • Ligand: ATP-binding; calcium; metal-binding; nucleotide-binding; zinc	1.42	Increase4, 8
Ephx2*	• Molecular function: hydrolase • Biological process: aromatic hydrocarbons catabolism; detoxification; lipid metabolism • Ligand: magnesium; metal-binding	− 0.71	Decrease4
Fabp7	• Biological process: transport • Ligand: lipid-binding	− 0.39	Decrease 2,5
Fkbp5	• Molecular function: chaperone; isomerase; rotamase	1.81	Increase 3,4,5,7,8
Gabra3	• Molecular function: chloride channel; ion channel; ligand-gated ion channel; receptor • Biological process: ion transport; transport • Ligand: chloride	− 0.81	Decrease6
Gfra1*	• Molecular function: receptor	1.26	Increase4, 8
Grm3*	• Molecular function: G-protein coupled receptor; receptor; transducer	− 0.63	Decrease4
Hpcal4*	• Ligand: calcium; metal-binding	− 0.69	Decrease4, 8
Hsph1	• Biological process: stress response • Ligand: ATP-binding; nucleotide-binding	− 0.80	Decrease3
Itm2a	N/A	− 0.73	Decrease 5
Mecp2	• Molecular function: DNA-binding; repressor • Biological process: transcription; transcription regulation	− 0.66	Decrease1, 4
Msmo1 / Sc4mol	• Molecular function: oxidoreductase • Biological process: lipid biosynthesis; lipid metabolism; steroid biosynthesis; steroid metabolism; steroid biosynthesis; sterol metabolism • Ligand: iron; NAD	− 0.77	Decrease 5
Myo1b*	• Molecular function: actin-binding; calmodulin-binding; motor protein; myosin • Ligand: ATP-binding; nucleotide-binding	1.32	Increase4, 8
Pcsk1	• Molecular function: hydrolase; protease; serine protease • Ligand: calcium	− 0.74	Decrease6
Pdia4*	• Molecular function: isomerase	− 0.71	Decrease4, 8
Plagl1	• Ligand: metal-binding; zinc	1.37	Increase 5
Prkcg*	• Molecular function: kinase; serine/threonine-protein kinase; transferase • Biological process: biological rhythms • Ligand: ATP-binding; calcium; metal-binding; nucleotide-binding; zinc	− 0.82	Decrease8
Pvalb	• Molecular function: muscle protein • Ligand: calcium; metal-binding	1.27	Increase 2
Pygm	• Molecular function: glycosyltransferase; transferase • Biological process: carbohydrate metabolism; glycogen metabolism • Ligand: nucleotide-binding; pyridoxal phosphate	1.29	Increase6
Qdpr*	• Molecular function: oxidoreductase • Biological process: tetrahydrobiopterin biosynthesis • Ligand: NADP	− 0.74	Decrease4
Rnpep*	• Molecular function: aminopeptidase; hydrolase; metalloprotease; protease • Ligand: metal-binding; zinc	− 0.78	Decrease4
Sgk1	• Molecular function: kinase; serine/threonine-protein kinase; transferase • Biological process: apoptosis; stress response • Ligand: ATP-binding; nucleotide-binding	1.37	Increase 2,3
Slc24a4*	• Biological process: antiport; calcium transport; ion transport; olfaction; potassium transport; sensory transduction; sodium transport; symport; transport • Ligand: calcium; potassium; sodium	1.58	Increase4
Slc9a3r1*	• Biological process: Wnt signaling pathway	− 0.74	Decrease8
Sun2*	• Biological process: meiosis	1.42	Increase8
Ugp2*	• Molecular function: nucleotidyltransferase; transferase • Ligand: magnesium; metal-binding	− 0.73	Decrease4
Zmat4*	• Molecular function: DNA-binding • Ligand: metal-binding; zinc	1.5	Increase8

First column represents the significant DE gene identified in this study, along with the UniProt function keyword(s) for each gene (second column) and the Log₂ fold change expression of the gene in our data set (third column). DE genes with an asterisk (*) represent genes identified in supplemental material from Chahrour et al. [4] and Veeraragavan et al. [34]. All function keywords were found from the UniProtKB/Swiss-Prot database [131] per respective DE gene. DE genes without a published function keyword in the UniProtKB/Swiss-Prot database are represented as “N/A.” The last column represents the previously identified RTT hits’ fold change expression direction, with each superscript representing that gene’s respective references. Superscript references are as follows: (1) Amir et al. [2], (2) Tudor et al. [20], (3) Nuber et al. [21], (4) Chahrour et al. [4], (5) Urdinguio et al. [24], (6) Ben-Shachar et al. [25], (7) Lin et al. [33], and (8) Veeraragavan et al. [34]. Rows in bold font represent previously identified RTT hits that were found to be specifically DE in the cortex

We then focused on comparisons of our RNA-Seq data to microarray studies also performed in Mecp2-deficient cortex [20, 24]. We specifically focused on these two studies for a number of reasons. First, studies have shown brain tissue region and cell type-specific gene and protein expression differences [30, 57–60]. While we do not negate the importance of other transcriptomic studies performed in the RTT field, we chose to compare our data set to studies performed in the same MeCP2 animal model, tissue region, and time point. A RNA-Seq study by Li et al. meets this criteria [9]; however, this study only reported upregulated genes and utilized different bioinformatics approaches. We found 175 upregulated genes in the Li et al. study were also present in our data set. Of these, 84/175 genes were also significantly upregulated in our data set (p < 0.05). We found 5 out of the 12 genes identified in the Tudor et al. microarray study and 6 out of the 20 genes identified in the Urdinguio et al. study were also significant DE genes in our data set (Fig. 1c). A common statistical significance threshold in RNA-sequencing studies is a q-value less than 0.05, which is meant to reduce the number of false positives observed, but may also underestimate biologically relevant changes. Reducing the significance threshold stringency to p < 0.05 in our data set, we identified an additional 2 genes from the Tudor et al. (7/12) and 2 genes from the Urdinguio et al. (8/20) study. Our RNA-Seq data correlates better with the individual microarray studies than these two studies correlate to each other, which is notable given strain differences between each of the individual studies, and that the Urdinguio et al. study utilized tissue pooled from three brain regions from the Mecp2 ^tm1.1Bird/y mouse model [24]. These genes and their putative functions as well as direction of expression levels for each hit are shown in Table 1. We identified an additional 19 gene targets previously found to be RTT gene hits in our data set that had a p-value of less than 0.05, but did not meet the false discovery rate, or FDR/q-value, cutoff of 0.05 (Additional file 2). Collectively, this suggests that our data is in concordance with previous literature on Mecp2-deficient mice.

In addition to known RTT hits, we also found a number of significant DE genes that have been previously identified as autism spectrum disorder (ASD) candidate genes (q < 0.05). For example, mutations in Auts2 [61, 62], Shank1–3 [63], and Foxp1 [64, 65] have been found in various ASD patient cases. In our data set, Auts2 and Foxp1 are significantly increased by 0.49 and 0.29 fold, respectively, while Shank1–3 are all decreased (Shank1: − 0.53 fold, Shank2: − 0.34 fold, Shank3: − 0.37 fold) (Table 1, Additional file 1). The gene expression changes observed in our data are robust, as many were identified across multiple models, different brain regions, and numerous species. Importantly, this also provides validation in utilizing this data set as an appropriate resource for comparing RNA-protein expression similarities in the RTT model.

Recent work demonstrated that MeCP2 selectively represses gene expression of long genes [30, 32]. These two studies also showed long genes are enriched in neurons and hypothesized that neuronal dysfunction in RTT could be in part due to repression of long gene expression [30, 32]. Recent studies have also indicated that disruption of long genes could also contribute to ASD pathogenesis [66–68]. For example, one study has shown that a heterozygous knockout of the autism candidate gene Chd8 led to an increase in dysregulated long gene expression [66], while another study found many ASD candidate genes are long genes [67]. Interestingly, recent studies have shown that topoisomerase inhibitors can regulate the expression of long genes in different ASD models [67, 68], suggesting that topoisomerase inhibitors could serve as a potential therapeutic for autism. Collectively, these studies highlight the need to understand the relationship between long gene expression and ASD and other ASD-related disorders such as RTT.

Given that MeCP2 is expressed in all CNS cell types, we extended earlier studies exploring the relationship between gene size and fold change expression within cell type-specific DE genes in our data. Accordingly, we next evaluated if the 391 identified significant DE genes in symptomatic Mecp2 ^Jae/y mice were enriched in a specific CNS cell type using a publicly available database of purified cortical CNS cell types [69]. The criteria for this comparison were as follows: (1) a gene which demonstrated a 3-fold or greater enrichment in 1 cell relative to all others was considered cell type enriched and (2) all other genes that varied less than 3-fold in any one cell type were considered to be non-specific. Using these criteria, 68% of the DE genes were not enriched in any cell type (Fig. 2a, Additional file 3). The other 32% of DE genes segregated into one enriched cell population. Of the DE genes associated with a specific cell type, 31 were enriched in neurons, 46 in astrocytes, 10 in microglia, 26 in oligodendrocytes (includes oligodendrocyte precursor cells, newly formed oligodendrocytes, and myelinating oligodendrocytes), and 14 in endothelial cells (Fig. 2a, Additional file 3). Additional comparisons to a second CNS cell type-specific database [60] further confirmed gene expression is disrupted across multiple cell populations in the symptomatic RTT brain (Additional file 4).

Based on the criteria outlined above, we cannot rule out the possibility that a gene we termed “not enriched in any cell type” is dysregulated in neurons. For example, the gene Camk1 has strong associations with neuronal transcription and synaptic activity. However, based on our and the Zhang et al. criteria [69], Camk1 is categorized as a “non-cell type enriched” gene. The non/cell type-specific hits may be a direct result of MeCP2 deficiency in neurons, MeCP2 deficiency in other CNS cell populations, or alternatively are downstream of MeCP2 binding in one or multiple CNS populations. We also cannot exclude the possibility that alterations observed in gene and protein expression in non-neuronal cell populations are due to the disease severity at the time point these experiments were performed. However, when we did a similar gene cell type-specific comparison between Chahrour et al.'s study and Zhang et al.'s database, we also observed a similar cell type-enriched distribution as the data presented in our study (data not shown), decreasing the possibility that our observed expression changes are due to disease severity.

We defined short genes as being less than 100 kb and long genes as larger than 100 kb based on prior studies [32] and calculated gene length based on the chromosomal location provided by the Cuffdiff differential gene expression output file. We found 81% of the total list of significant DE genes was short in our data (Fig. 2b, Additional file 5). Of these short genes, 237 had decreased and 79 had increased gene expression. Similar to previous work [30, 32], we found 71% (53/75) of the long genes exhibit increased gene expression (Fig. 2b, Additional file 5). Although the majority of genes in our data set are short and repressed, the neuronal-enriched DE genes in our data support previous studies showing a bias for over-expressed long genes in neuronal populations [30] (Fig. 2b, Additional file 5). Additionally, increased long gene expression was observed in oligodendrocyte-enriched and non-cell type-specific DE genes (Fig. 2b, Additional file 5). In contrast, the majority of astrocyte-enriched DE genes showed decreased gene expression, in which 92% (36/39) of short genes and all 7 long genes exhibited decreased gene expression (Fig. 2b, Additional file 5). Similar trends were observed in endothelial-enriched DE genes, where 71% (10/14) of the genes were short and decreased (Fig. 2b, Additional file 5). Furthermore, within the microglial-enriched DE gene group, all 10 genes are short and have decreased gene expression (Fig. 2b, Additional file 5). Together, these data provide an additional dimension of cell type-specific gene expression changes in Mecp2 ^Jae/y animals.

Protein expression in symptomatic Mecp2 ^Jae/y whole cortex

To gain a more complete view of cellular and molecular dysfunction in RTT, we examined global protein abundance changes in symptomatic Mecp2 ^Jae/y whole cortex lysates using a DIA LC-MS/MS proteomics approach (see the “Methods” section; [50, 53]). We identified 465 significant protein changes out of 4789 total quantified proteins. Of the 465 significant proteins, 299 proteins had increased fold changes in RTT mice while 166 had decreased fold changes compared to WT (Fig. 3a, Additional file 6). To demonstrate the robustness and sensitivity of the expression changes observed in our proteomics data set, we performed Western blot analysis on 3 top-ranking proteins based on large fold change differences (fold change greater than ± 3; MeCP2, CIRBP, RNPEP) along with 4 proteins with smaller fold change differences (fold change less than ± 3; HSPH1, NEFM, QDPR, mGluR3) (Additional files 6 and 7). As expected, the MeCP2 protein was one of the most highly decreased and significant proteins in the data set and was also significantly decreased by Western blot analysis (Additional files 6 and 7). Western blot analysis also confirmed that the remaining 6 targets were significantly different in the Mecp2-deficient cortex (Additional file 7). Additionally, representative mass spectrometry chromatograms confirmed qualitative changes in peptides detected in WT versus Mecp2-deficient cortical lysates (Additional file 7). Collectively, these results lend support to the robustness and sensitivity of our proteomics data set.

We next examined whether any of the significant proteins were previously identified RTT hits using the same criteria as outlined in the RNA-Seq analysis section. Of the 465 significant proteins, 37 were identified as previously published RTT targets, with the majority identified from previous transcriptome studies (Fig. 3b, Additional file 8). Three of the 37 hits (clusterin, Apo-J/Clu; alpha-1-antitrypsin 1–1, AAT/Serpina1a; and apolipoprotein A-I, Apo-AI/Apoa1; Fig. 3b, Additional file 8) were identified from the limited proteomics studies published in RTT [36–39] and have decreased expression. Two significantly increased proteins in our data set (transcriptional activator protein Pur-alpha, Purα/Pura, ↑1.19 fold; and core histone macro-H2A.2, mH2A2/H2afy2, ↑1.43 fold; Fig. 3b, Additional file 8) were previously identified as MeCP2 interacting proteins [4] as was FK506-binding protein 51 (FKBP5; Additional file 8) [21].

Our proteome data indicate cellular dysfunction in RTT may be facilitated by disrupted proteostasis (Table 2). We observed significantly decreased abundance of multiple heat shock proteins (HSPH1, HSP90α, HSP90β, and HSPA4), including the small heat shock proteins alpha-crystallin B chain (CRYAB) and HspB1 (Table 2). In support of these findings, a recent study has shown that microglia isolated from pre-symptomatic RTT female mice show decreased expression in heat shock genes [73]. Because of the important role that heat shock proteins play in regulating protein folding and stability, it is possible that expression changes in these heat shock proteins could be influencing the protein expression changes we observe in our proteomics data set. Furthermore, the protein IMPACT, which mediates translation in response to numerous cell stressors, was also upregulated in the proteome data set (1.8-fold; Table 2). Elevated physiological stress can also affect the protein expression levels of FKBP5 (elevated in response to glucocorticoid signaling [21]), the diacylglycerol-binding protein DGKγ, and the guanine nucleotide exchange factor RasGRF1 [74, 75], all of which are significantly increased in our data (Additional file 6).

Recent reports indicate altered metabolism in RTT. We observed significant expression of proteins involved in ketone body metabolism and utilization (AACS) as well as fatty acid oxidation (ACADVL) and glycogen metabolism (UGP2) (Table 3). Lanosterol synthase (LSS), a key enzyme in cholesterol synthesis, was elevated 8.1-fold in our data set (Table 3). Cholesterol synthesis has been shown to be disrupted in a suppressor screen in Mecp2-deficient mice [76]. It is conceivable that elevated LSS protein expression contributes to elevated cholesterol triglycerides and LDLs commonly observed in RTT patients [77, 78]. Our data also indicate altered monoamine metabolism as evidenced by significant downregulation of two key enzymes (QDPR, SPR) involved in the production of dopamine, norepinephrine, and serotonin (Table 3). The product of these enzymes is tetrahydrobiopterin (BH4), an essential cofactor in the production of amine neurotransmitters. Furthermore, amine oxidase A (MAOA), which catalyzes the oxidative deamination or degradation of these same neurotransmitters, is upregulated in the proteome data set (Table 3). Together, these findings provide additional support for prior studies indicating that dopamine, norepinephrine, and serotonin are downregulated in RTT mice and patients [79, 80].

Finally, we identified a number of significantly altered proteins associated with S-adenosylmethionine-dependent methylation (GAMT, SPDSY, and MTHFD1L; Table 3). This particular pathway converts methionine to S-adenosylmethionine (SAM), whereby the methyl group from SAM is transferred to different substrates such as DNA [81, 82]. This finding is of interest to RTT because MeCP2 binds methylated and un-methylated genomic regions to regulate gene transcription [4], and global MeCP2 methylation-binding patterns in RTT have been well documented [13, 31, 32, 83, 84]. An additional methyltransferase also had significant differential abundance between RTT and WT animals (PRMT5; Additional file 6). Methyl donor groups that are utilized by SAM-dependent methylation mechanisms include creatine, folate, folinic acid, and betaine, all of which have also been implicated in RTT [85]. Collectively, this suggests aberrant regulation of methyl donor group availability to properly regulate DNA methylation, and consequent downstream gene expression, in RTT.

Significant proteins in symptomatic Mecp2 ^Jae/y whole cortex are expressed in multiple CNS cell types

To determine if any of the significantly abundant proteins identified in our proteomics data set were also enriched in a specific cell type, we compared the 465 significant proteins to the same CNS cell type database as the RNA-Seq data. Similar to our DE gene findings, most significantly altered proteins were non-cell type specific (Fig. 3c, Additional file 3). Of the cell type-specific significantly regulated proteins, 12 were neuronal-enriched, 14 astrocyte-enriched, 4 microglial-enriched, 34 oligodendrocyte-enriched, and 10 endothelial-enriched (Fig. 3c). Similar to the RNA-Seq data, the proteomics data also suggest that loss of MeCP2 has downstream effects on protein expression in different CNS cell populations.

It should be noted that in our study, we have solely examined gene/protein expression in whole cortical tissue homogenate. A number of studies have identified this “dilution effect” [30, 86, 87], where certain genes in one cellular population could be masked by the gene expression from a heterogeneous cellular population. Additionally, gene and proteins typically thought of as “cell type specific” may lose their specificity in the RTT brain. Future studies examining gene and protein expression from multiple cell types simultaneously isolated from the same RTT animal model and tissue region are needed to properly address these potential confounds. We also cannot exclude the possibility that changes in gene and protein expression could be due to differences in cell population number. Future studies utilizing quantitative stereology and multi-channel immunohistochemistry should be done to examine cell counts in WT versus Mecp2-deficient brain tissue.

Transcriptome-proteome expression comparison identifies novel hits

We next compared the proteome to the RNA-Seq data set. Previous studies making similar comparisons across species reported low correlation between these two types of data [88–91]. Our comparison across the transcriptomic and proteomics data resulted in low correlation (Pearson’s R of 0.12, Fig. 4a). Within the detected gene-protein matches, which we define as a detected gene having a corresponding detected protein, we identified 35 significant gene-protein matches (gene expression q < 0.05, protein expression p < 0.1) with a Pearson’s R of 0.74 (Fig. 4b and Table 4). When we relaxed the stringency in the DE gene list to p < 0.05, this value increases to 77 significant gene-protein matches with a Pearson’s R of 0.67 (Additional file 9).

Gene/protein	UniProt accession	Gene fold change	Protein fold change	RTT hit?	Putative function
Mecp2/MeCP2	Q9Z2D6	− 0.66	− 10.6	Yes1, 2	DNA binding protein; transcriptional regulator
Rnpep/RNPEP*	Q8VCT3	− 0.78	− 3.47	Yes2	Peptide catabolic process
Pdia4/ERp72*	P08003	− 0.71	− 2.22	Yes2,3	ER protein processing
Hsph1/HSPH1	Q61699	− 0.80	− 2.13	Yes4	Protein chaperone
Rnpep/RNPEP(Ac)*	Q8VCT3	− 0.78	− 2.02	Yes2	Peptide catabolic process
Aldh1a1/ALDH1A1*	P24549	− 0.67	− 1.95	Yes2	Retinal dehydrogenase
Slc9a3r1/NHERF1(P)*	P70441	− 0.74	− 1.79	Yes3	Scaffold protein
Ugt8a/UGT8	Q64676	− 0.79	− 1.76	---	Sphingolipid metabolism
Qdpr/QDPR (or DHPR)*	Q8BVI4	− 0.74	− 1.68	Yes2	Monoamine metabolism (tetrahydrobiopterin biosynthesis)
Ugp2/UGP2*	Q91ZJ5	− 0.73	− 1.67	Yes2	Glycogen synthesis
Car2/Ca2	P00920	− 0.72	− 1.61	---	Reversible hydration of carbon dioxide
Grm3/mGluR3*	Q9QYS2	− 0.63	− 1.6	Yes2	Protein coupled glutamate receptor
Ephx2/EPHX2*	P34914	− 0.71	− 1.59	Yes2	Cholesterol synthesis
Hpcal4/HPCAL4*	Q8BGZ1	− 0.69	− 1.51	Yes2,3	Calcium binding
Spr/SPR	Q64105	− 0.77	− 1.47	---	Monoamine metabolism (tetrahydrobiopterin metabolic process)
Sash1/SASH1	P59808	− 0.82	− 1.44	---	p38 MAPK and NIK/NF-kappaB signaling
Calr/CRT*	P14211	− 0.80	− 1.43	Yes3	ER calcium-binding protein
Aacs/AACS*	Q9D2R0	− 0.79	− 1.37	Yes2	Fatty acid and lipid metabolism
S1pr1/S1P1	O08530	− 0.78	− 1.28	---	Cell-cell adhesion
Slc9a3r1/NHERF1(Ac)*	P70441	− 0.74	− 1.24	Yes3	Scaffold protein
Nova1/NOVA1	Q9JKN6	1.26	1.28	---	RNA/mRNA-binding protein
Slc24a4/Nckx4*	Q8CGQ8	1.58	1.3	Yes3	K+-dependent Na+/Ca2+ exchanger
Kcnab3/Kvβ3.1	P97382	1.27	1.43	---	Voltage gated K+ channel
Me3/NADP-ME3	Q8BMF3	1.28	1.46	---	Mitochondrial NADP(+)-dependent malic enzyme
Tfrc/TfR	Q62351	1.32	1.48	---	Cell surface receptor required for iron uptake
Dgkg/DGKγ*	Q91WG7	1.42	1.51	Yes2,3	Lipid metabolism
Prkcg/PKCγ(P)*	P63318	− 0.82	1.66	Yes3	Protein kinase; LTP
Rasgrf1/RasGRF1	P27671	1.32	1.68	---	Stimulates the dissociation of GDP from RAS protein
Zmat4/ZMAT4*	Q8BZ94	1.50	1.88	Yes3	DNA/RNA and p53 binding
Myo1b/Myo1b*	P46735	1.32	2.02	Yes2,3	Actin binding; calmodulin binding
Fkbp5/FKBP5	Q64378	1.81	2.02	Yes2–6	Immunoregulation; protein trafficking
Wipf3/WIPF3	P0C7L0	1.23	2.26	---	Cytoskeleton organization
Itih3/ITI-HC3	Q61704	1.52	2.37	---	Anchor protein between hyaluronan and matrix proteins
Gfra1/GDNFRα1*	P97785	1.26	2.64	Yes2,3	RET and RAF/MAP kinase signaling cascade
Sun2/SUN2*	Q8BJS4	1.42	2.86	Yes3	Nuclear envelope protein

Significant genes were defined as having a q < 0.05, and significant proteins having a p < 0.1. Gene and protein names are provided along with the UniProt accession number, gene and protein fold change expression values (respectively), whether the gene-protein match has been previously identified as a RTT hit, and the respective putative function according to the UniProtKB/Swiss-Prot database [131]. DE gene-protein matches with an asterisk (*) represent genes/proteins identified in supplemental material from Chahrour et al. [4] and Veeraragavan et al. [34]. Gene-protein matches with “(Ac)” denote that the corresponding protein had a significant acetylation PTM, and a “(P)” denotes a significant phosphorylation PTM. For the RTT hit column, a “---” indicates the match has not been identified as a RTT hit prior to this study. Superscripts in the RTT hit column denote references for matches identified as a RTT hit and are as follows: (1) Amir et al. [2], (2) Chahrour et al. [4], (3) Veeraragavan et al. [34], (4) Nuber et al. [21], (5) Urdinguio et al. [24], and (6) Lin et al. [33]

We found 23 of the 35 significantly DE gene-protein matches were previously identified as RTT RNA regulatory targets (Table 4). The remaining 12 matches have not been described in previous RNA studies (Table 4). Within the gene-protein matches with similar regulatory patterns, Ephx2/EPHX2, Dgkg/DGKγ, Me3/NADP-ME3, Qdpr/QDPR, Slc24a4/SLC24A4, Ugt8a/UGT8, Aacs/AACS, Rnpep/RNPEP, and Spr/SPR are implicated in metabolic pathways (Table 4) [92–100]. We also found that Itih3/ITI-HC3, Slc9a3r1/NHERF1 (both phosphorylated and acetylated modifications), and Rnpep/RNPEP are associated with protein scaffolding/stability (Table 4) [99, 101–103]. Additionally, S1pr1/S1P1, Calr/CALR, Hpcal4/HPCAL4, and Rasgrf1/RasGRF1 are implicated in calcium-mediated processes (Table 4) [104–107].

Pathway analysis in transcriptome and proteome implicate similar biological pathways

In an attempt to understand the underlying biology of disrupted genes and proteins, we utilized Ingenuity® Pathway Analysis (IPA) on the protein and RNA data sets. For a comprehensive list of all identified pathways in the individual transcriptome and proteome data sets as well as the corresponding genes and proteins in each pathway, see Additional file 10. To exemplify the diversity of identified significant pathways, we highlight functional categories that are similar or shared between the transcriptome and proteome data sets in Figs. 5 and 6. First, we identified pathways associated with general cellular and molecular dysfunction. This included categories such as cell cycle, cell components/structure/general function, and lipids and metabolism (Fig. 5a-c, Additional file 10). For the cell cycle category, 36% (4/11) of the identified pathways are associated with G₂ and G₂/M phase, which were found exclusively in the proteomics data set (Fig. 5a, Additional file 10). Over 80% of the proteins associated with these specific pathways had significantly increased expression, suggesting hypersensitivity and/or an accumulation of cells in the G₂/M phases, a phenomenon previously observed in RTT primary human fibroblasts [108–110].

Within the cell components/structure/general function category, we found a number of pathways related to cell death, tissue/cellular morphology, and G-protein signaling shared between the transcriptome and proteome (Fig. 5b, Additional file 10). In the RTT brain, the primary morphologic change in humans is reduced brain size [111]. Additional morphological changes include simplified dendritic arborizations and altered spine density and morphology [41, 112]. A recent study also identified reduced astrocyte branching and overall process length [72]. Supporting this, pathway analysis identified several disrupted pathways associated with tissue and cellular morphology such as “morphology of cerebral cortex”, “abnormal morphology of cells”, “formation of intermediate filaments”, “quantity of plasma membrane projections”, and “differentiation of cells” (Fig. 5b, Additional file 10). Additionally, pathway analysis supports previous findings that metabolites including glucose [113–115], lipids/phospholipids [77, 78, 80, 115–117], TCA cycle intermediates [115], biogenic amines [79, 80], glutamate [80, 117], creatine [80, 85], and steroids such as cholesterol and glucocorticoids [76, 78, 114, 118–120] are aberrantly regulated in RTT (Fig. 5c, Additional file 10).

Recent interest in different neural cell types in RTT prompted us to examine dysfunction in cell type-specific pathways. Given the wealth of information on neuronal dysfunction in RTT, we unsurprisingly identified disrupted pathways associated with synaptic/neurotransmission, neuronal morphology, and development. Proteome pathway analysis also indicated disrupted conduction and neuronal structure/organization (Fig. 6a). These pathways and expression trends are consistent with previously published work in RTT [41, 111]. Within the glial functions group, we found the majority of pathways in both the transcriptome and proteome data sets are associated with myelination processes, along with astrocyte morphology, reactivation, and apoptosis (Fig. 6b). A high number of oligodendrocyte function-based pathways were observed in the proteome (8/9 pathways) in which 7 of the pathways have decreased protein expression (Fig. 6b). While the exact mechanisms of how oligodendrocytes contribute to RTT pathogenesis have received little attention, previous studies have implicated an abnormal oligodendrocyte/myelin involvement in the disease [17, 121–123]. Altogether, the glial pathways seem to suggest abnormal glial morphology as well as aberrant myelination functions in RTT.

We additionally identified unique pathways associated with inflammation/immunology and blood vessel/vasculature pathways as disrupted in the transcriptome and proteome data sets (Fig. 6c, d, respectively). Within the inflammation/immunology category, 90% (9/10) of the pathways had decreased gene and protein expression (Fig. 6c). Furthermore, 7 immunology/inflammation pathways were found to be shared in common with both the transcriptome and proteome data sets (Fig. 6c, Additional file 10). Indeed, our pathway analysis has identified pathways such as “T cell development”, “abnormal quantity of cytokine”, “recruitment of lymphocytes”, and “activation of lymphocytes” (Fig. 6c), in which disruptions in T cell/T lymphocyte differentiation and cytokine regulation have been documented in RTT [124, 125]. These identified immunology/inflammation pathways highlight a disrupted inflammatory component in RTT, a topic that has received considerable recent attention [124, 125].

Finally, we observed a number of pathways associated with blood cells such as “production of blood cells”, “expansion of blood cells”, “recruitment of blood cells”, and “elimination of blood cells” (Fig. 6d). This is particularly interesting as platelet defects have been observed in RTT patients with FOXG1 mutations [126]. While few studies have examined general vasculature in RTT, defects in vascular function have been previously observed [127]. The above identified pathways suggest a dysregulated function of the vasculature, and perhaps an increased risk of blood-brain barrier breakdown (a pathway also identified in the proteomics data, Fig. 6d), in RTT symptomatic WCX.

Why the proteome?

Much of our understanding of MeCP2 function is derived from transcriptomic studies, with the general assumption that alterations in the transcriptome correlate with proteomic changes. Challenging this central dogma, recent studies indicate that RNA-protein expression correlations are low even when comparing the same tissue samples [88–91]. There are a number of reasons given for the observed discrepancies, including differences in stability and lifetime of the two types of molecules, posttranscriptional and posttranslational modifications, and protein turnover [90, 91]. These changes are likely to be exacerbated in a disease state, particularly a disease such as RTT where data indicates RNA metabolism and proteostasis are disrupted (Table 2). Thus, the primary benefit of performing a direct comparison between transcriptomic and proteomic data is the ability to complement the knowledge gained from each of these omics technologies, providing a more holistic understanding of the interplay between gene and protein expression regulation. The benefits of such comparisons are reflected in our pathway analyses, where “neuronal transmission”, “synaptic transmission”, “excitation of neurons”, “morphology of neurons”, and “development of neurons” pathways are unique to the transcriptome data set. Supplementing these transcriptomic pathways, pathways identified using the proteome data indicated “abnormal conduction”, “damage of synapses”, “atrophy of axons”, and “organization of neurofilaments” as disrupted. These findings indicate that evaluating data obtained from multiple omics platforms provide a more complete understanding of disease biology.

While the current study has only examined transcriptomic and proteomic expression changes in symptomatic Mecp2-deficient male mice, we recognize that utilization of a heterozygous female Mecp2-deficient mouse model is more clinically relevant. To date, the majority of transcriptomic and proteomic studies have been performed in male animal models. Therefore, we chose to use symptomatic Mecp2-deficient male mice for comparison with existing data in the RTT field. Future studies should be directed toward understanding proteomic and transcriptomic changes occurring in the same RTT female mouse model and tissue region.

Cell type-specific gene and protein expression