AI-enabled drug prediction and gene network analysis reveal therapeutic use of vorinostat for Rett Syndrome in preclinical models

Rett Syndrome (Rett) is a neurodevelopmental disorder that is also known to clinically impact multiple organs and systems in the body. It is primarily caused by mutations of the gene encoding methyl-CpG binding protein (MeCP2), a ubiquitous nuclear protein that acts as a transcriptional repressor and activator affecting hundreds of genes across the genome¹. Rett is characterized by a high degree of cellular heterogeneity, with a mixture of wild type and MeCP2 mutant cells within the tissues of the same individual² that leads to subtle but widespread gene misregulation³. Clinical criteria and research on Rett have primarily focused on neuromotor, social, and cognitive impairments and other CNS symptoms (e.g., seizures, repetitive purposeless movements, regression, microcephaly, loss of speech, autistic features)^4,5. However, these patients also commonly exhibit multiple clinical symptoms outside of the CNS, including respiratory dysfunction, gastrointestinal issues, systemic inflammation, disrupted metabolism, circadian rhythm and sleep disturbances, abnormal bone density, kidney function, nociception derangement, which are a major source of morbidity that contribute to a shortened lifespan^6,7. It is important to note that although Rett is considered a neurodevelopmental disorder and therapeutics development has focused on CNS-related endpoints, these system-wide symptoms are often the most disturbing to patients and families, and thus there is a need to ameliorate symptoms of this disease across the entire body.

Monogenic diseases are theoretically amenable to gene therapy treatment, and Rett has been shown to be reversible in a mouse model of X chromosome reactivation⁸ and post-natal expression of MeCP2^9,10. However, even if current challenges around delivery efficacy and safety could be overcome, Rett and other X-linked disorders require gene dosage compensation to avoid toxicity due to overexpression of the corrected gene because even a modest 20% deviation of MeCP2 expression levels can lead to Rett symptoms^11,12. Indeed, several clinical trials have been initiated to test gene therapies that regulate MeCP2 gene expression (NCT05898620, NCT06152237, NCT05606614). However, due to the limited number of patients enrolled and uncertain outcomes, effective pharmaceutical treatments for Rett and other neurodevelopmental disorders are still needed. Trofinetide, a synthetic IGF-1 C-terminal tripeptide activating BDNF expression, was approved by the FDA to treat Rett syndrome in March 2023, but it has only shown moderate efficacy and its use is limited by significant gastrointestinal adverse events¹³. Sarizotan, a 5-HT1A and D2 receptor agonist with promising efficacy in mouse models of Rett¹⁴, was recently terminated after clinical trial results failed to demonstrate efficacy (NCT02790034). A Phase III study with blarcamesine also failed to meet its endpoints and was terminated (NCT03941444). Esketamine, an N-methyl D-aspartate receptor (NMDAR) agonist, is also being evaluated in the clinic (NCT03633058) due to its reduction of neuroexcitatory glutamate. But these and other therapeutic paths for Rett pursued thus far target only a subset of patients’ CNS symptoms and have not shown any effects on extra-CNS morbidities.

Here we combine an integrated drug repurposing approach using a computational model that combines artificial intelligence (AI) and gene expression network analysis to predict existing drugs that may be used for treatments¹⁵ with phenotypic screening in an phenotypically diverse in vivo model of Rett syndrome generated using CRISPR technology in Xenopus tadpoles to assess whole-body efficacy for clinically-relevant metrics. Our results identify vorinostat as a lead candidate, and we demonstrate that it is highly effective in reversing behavioral changes and multiple organ abnormalities in a mouse Rett syndrome model even when treatment was initiated after the onset of symptoms. These findings demonstrate the utility of this target-agnostic platform for predicting drugs to restore multi-organ system function in Rett syndrome and potentially other complex diseases, in addition to leading to a potential therapeutic option for this devastating neurological disease.

X. laevis care and rapid generation of Rett syndrome tadpole models

All animal experiments and procedures were reviewed and approved by the Harvard Medical School (HMS) Institutional Animal Care and Use Committee regulations (Protocol #IS2714). Xenopus laevis embryos and tadpoles were housed at 18 °C with a 12/12 h light/dark cycle in 0.1X Marc’s Modified Ringer’s (MMR) medium. Tadpoles were tested up to Stage NF48 and thus are sexually undifferentiated, since the sex-determining period starts on Stage NF50¹⁶. Tadpoles were fed 3 times per week with Sera Micron Nature fry food. Xenbase¹⁷ (www.xenbase.org) was used to search for MeCP2 homolog in Xenopus laevis (Xenopus laevis J-strain 9.2 version). MeCP2.L (XB-GENE-17346751), MeCP2.S (XB-GENE-494744), found on Chromosome 8 L and 8S, show 75% and 70% sequence identity with human MeCP2, respectively. An important biological distinction is that while human MeCP2 is found on the X chromosome, Chromosome 8 is not part of the sex determination alleles of homomorphic chromosome pairs in X. laevis. Cas9 target sites were selected using CHOPCHOP¹⁸ on the X. laevis J-strain 9.2 version, which facilitated the selection of guide RNA with no predicted off-targets. We used the default settings to generate a list of target sequences. We selected 4 targets from this list for each of the L and S forms. This selection was based primarily on the ranking provided by CHOPCHOP and as a function of their ___location on MeCP2 to target exons coding for the methyl-CpG-binding ___domain (MBD, exons 2 and 3) and Transcriptional repression ___domain (TRD, exon 3). The 8 sgRNA sequences are presented in Supplementary Table 1 and were synthesized as modified sgRNA by Synthego Inc.

The sgRNA were resuspended to 100 μM in 0.1X Tris EDTA (pH 8.0). Then, an equimolar sgRNA mix was made using 2 μl of each CRISPR sgRNA, for a total of 16 ul. The Cas9 RNP complex was formed by mixing 75 pmol of the sgRNA mix with 75 pmol of Cas9 in annealing buffer (5 mM HEPES, 50 mM KCl, pH 7.5), for a total volume of 100 μL, then incubated at 37 °C for 10 min. The RNP was kept frozen at −20 °C until the day of the injection. To generate MeCP2 knockdown models of Rett syndrome, Xenopus embryos were fertilized and maintained at 14 °C until the 4-cell stage. Each cell was injected with ~2 nL Cas9 ribonucleoprotein complexed with sgRNA (Cas 9 RNP) per injection resulting in a final amount of 1.5 fmol of Cas9 RNP per injection in each of the 4 cells, and we did not see any adverse effect when controlling for Cas9 or sgRNA at that amount. Cas9 RNP contained an equimolar mix targeting 6 sites of the X. laevis MeCP2 gene on the chromosomes 8 L and 8S. Note that the 1.5 fmol of Cas 9 RNP used per injection is two to five time less than what has been commonly used in previously published studies for Xenopus¹⁹ as we found that higher amounts of sgRNA had adverse effects (not shown). Following injections, embryos were kept in 0.1X MMR at 18 °C as described above, for the 18 days duration of the experiment. This approach generated mosaic knockout models of MeCP2 knockdown, which presented with a heterogeneous degree of phenotypes and severities when assessed at the population level. This approach is versatile, such that increasing the concentration of Cas9 RNP or injecting at the 8-cell stage leads to a narrower range of heterogeneity, with mostly severely affected tadpoles and fewer viable embryos, reflecting the potential lethality of complete MeCP2 knockouts.

PCR and fragment analysis to validate MeCP2 gene disruption

MeCP2 editing efficiency was measured using Indel Detection by Amplicon Analysis (IDAA)²⁰, which leverages PCR, using a pair of primers and universal fluorescein labeled oligo (Supplementary Table 2).The primers were synthesized by IDT. The amplicons were purified according to the manufacturer’s instructions (Zymo Research) and submitted for fragment analysis using capillary electrophoresis on an ABI 3730XL (performed by Genewiz), where each sample ran with an internal DNA standard. ProfileIT IDAA analysis software (COBO Technologies) was used to analyze the fragments and visualize the distribution of indels and visualize out-of-frame indels²¹. IDAA is a robust method with 1-bp resolution suited for the analysis of a tetraploid genome without the need for a deep sequencing approach.

We analyzed the indels of 3 target sites on MeCP2.S of 14 tadpoles presenting with seizures and 3 non-edited tadpoles. Using a window of ±100 bp from the target site, we observed varied editing efficiencies depending on the target site. The likely out-of-frame indels ratio averaged 47% for Xl-mecp2s-g01 on exon 3, with the lowest being 12% for a given tadpole and highest being 95%, 48% for Xl-mecp2s-g02 on exon 2 (min/max 20–84%), and ~82% for Xl-mecp2s-g04 on exon 2 (min/max 7–91%). Overall, 98% of indels where between 1 and 30 bp, and 84% were deletions.

Microarray analysis

Tadpoles were sacrificed at stage 50 of development and lysed using a QIAGEN TissueLyser II bead mill (30 Hz, 2 × 30 s) in a 2 mL tube with a 2.4 mm steel beads (Omni) and 1 mL phosphate buffered saline buffer at 4 °C. 250 mL of lysate was used for RNA extraction using a QIAGEN RNeasy mini kit. RNA samples were DNAse treated. Xenopus laevis Genome 2.0 Array (Affymetrix) microarrays were processed by Advanced Biomedical Laboratories (Cinnaminson, NJ). RNA samples were processed using a Nugen Ovation PICO WTA System V2 kit. The resulting cDNAs were purified using a Qiagen MinElute PCR Purification Kit following the modifications outlined in the Nugen protocol. The cDNAs were fragmented and labeled using a Nugen Encore Biotin Module. Hybridization solutions were prepared by combining the fragmented, biotin-labeled cDNAs with hybridization cocktail (Affymetrix Hybridization, Wash, and Stain Kit). The mixtures were incubated in a thermal cycler at 99 °C for 2 min followed by 45 °C for 5 min. then loaded on Xenopus laevis Genome 2.0 arrays and incubated for 16–20 h at 45 °C and 60 rpm in an Affymetrix Hybridization Oven 645. Following hybridization, arrays were washed and stained on Affymetrix Fluidics Station 450 s using the Affymetrix FS450_0001 protocol with the stains and buffers supplied in the Affymetrix Hybridization, Wash, and Stain Kit. The stained arrays were scanned at 532 nm using an Affymetrix GeneChip Scanner 3000.

Transcriptomics data processing

Microarray datasets were RMA normalized. RNAseq gene counts were normalized using size factors similar to the normalization performed by the R package DESeq. Genes missing values were removed and duplicates were combined by taking the maximum across counts. Gene expression heatmaps and hierarchical clustering plots were generated for all groups using the R package limma (version 3.42.2) and Python package seaborn (version 0.9.0). Volcano plots were generated using the Python package bioinfokit (version 0.7).

nemoCAD drug prediction algorithm

In the present study, we used the nemoCAD computational tool that combines AI and gene network analysis to predict drugs that would shift MeCP2-edited tadpoles to a control tadpole state based on transcriptome state for each condition. This algorithm was previously developed and use to identify existing drugs that could be repurposed for treatment of COVID-19 where it correctly predicted the impact of only a subset of statin compounds on patient survival¹⁵. The nemoCAD tool utilizes machine learning and pre-computed interaction probabilities of drug-gene and gene-gene interactions based on analysis of differential gene expression signatures of a disease state and appropriate controls to identify compounds capable of changing a transcriptional signature indicative of one biological state to another state of interest (e.g., reverting a disease state to a healthy state).

The NemoCAD repurposing tool is first used to identify transcriptome-wide differential expression profiles between two biological states in the input transcriptomic dataset (experimental or published) and to define the target normalization signature, i.e., the subset of genes whose expression levels need to be reversed in order to revert one state to the other (e.g., normalize the diseased state). Pairwise analysis, that implicitly assumes genes are expressed independently of one another, is carried out based on the comparison between gene expression profiles of >19,800 compounds found in the LINCS database release v1²² and the target normalization signature. Multiple correlation statistics (e.g., Pearson correlation, cross-entropy) across all differentially expressed genes are calculated and a combination score (e.g., Pearson correlation divided by the cross-entropy) is computed.

Putative predictions that incorporate gene-gene dependencies are then generated separately based on Bayesian network analysis on a regulatory and drug-gene interaction network architecture defined using publicly available databases of gene-gene interactions based on single gene knockout datasets in human cells (KEGG, TRRUST) and reference transcriptional signatures of drugs (LINCS, CTD), as described below. By combining a directed unweighted network structure with interaction probabilities for the connecting network edges, the constructed network is a weighted directed graph consisting of all possible paths that connect at least 2 genes of interest from the relevant genes within the target transcriptomic normalization set and the drugs; it therefore encodes the entire region of influence of a given list of genes and the drugs that can reverse their gene expression profiles in the desired manner. This network is then used as input for a message-passing algorithm (e.g., loopy belief propagation algorithm^23,24) in which the marginal probability distributions of drugs being turned on given the expression state of every gene are computed using the joint probability distribution encoded in the marginal distributions. This network-based approach leads to a ranking of all the drugs based on the probability of them inducing the desired transcriptomic signature across all genes in the network, which is separate from the initial ranking based on correlation analysis described above. This method has the advantage of incorporating gene-gene dependencies, but its success relies on the quality of network structures curated in the source databases.

Ultimately, drugs are prioritized that have both high correlation with the desired transcriptional signature change, representing desired drug-gene interactions, and high probabilities of being able to reverse the diseased state predicted by the Bayesian network analysis, representing gene network-level drug effects that can counter network-level disease processes. nemoCAD also enables visualization of the architecture of the network and subnetworks, and thus provides insight into potential molecular targets that can be verified using additional experiments. nemoCAD also can be used to optimize treatment regiments by iteratively inputting animal omics data following drug treatment to further tune drug predictions and develop drug combinations that might synergize in unpredictable ways. The approach underlying nemoCAD avoids the needs for a priori target or pathway inputs and, in fact, enables their discovery. Trofinetide was not in the digital molecule database, so we were unable to make comparisons with our top drug candidates.

We used multiple methods to assess the robustness of the nemoCAD algorithm (Supplementary Fig. S1). To demonstrate resilience of the drug predictions with a range of potential input transcriptome signatures, microarray data were artificially degraded through the addition of stochastic noise while drug rank was read out over 20 independent simulations per noise level. Additionally, drug prediction results were benchmarked against CMAP predictions that rely on the same underlying LINCS drug signature datasets but lack the network analysis of nemoCAD. Lastly, confirmation of self-consistency was performed using the vorinostat-treated cell line expression data from LINCS and reading out the top drug predictions from nemoCAD.

Transcriptomics gene network construction

Gene regulatory networks were inferred from gene expression data using the Bioconductor package GEne Network Inference with Ensemble of trees (GENIE3) release 3.12 in RStudio (R version 3.6.2)²⁵. GENIE3 infers a weighted, directed gene regulatory network from expression data by decomposing the prediction of a network between k genes into k different regression problems. For each individual regression analysis, the expression pattern of one of the genes (target gene) is predicted from the expression patterns of all the other genes (regulator genes) using ensembles of regression trees. GENIE3 produces an adjacency matrix representation of the network with k nodes in which each node represents a gene, and an edge directed from one gene i to another gene j indicates that gene i regulates the expression of gene j. Based on known mechanisms involved in Rett Syndrome (Supplementary Table 3), a subset of genes was selected as candidate regulators: BDNF, FCRL2, FMR1, MeCP2, NTRK2, PER1, PER2, and PUM1. These potential regulators were calculated over all available target genes in the transcriptomics datasets (10,934 genes) for MeCP2 knock down and cas9-injected control samples. Digraph objects were constructed from adjacency matrices and visualized in MATLAB R2020a (Mathworks; Natick, MA) for each of the MeCP2 knock down and control groups. The strongest regulator gene-target gene relationships were highlighted by filtering each network at a consistent edge weight (e_w) threshold (e_w = 0.3) prior to plotting.

Cross-species transcriptomics analysis

Gene networks involved in Rett Syndrome in human patients and MeCP2 knock down in laboratory animal models were compared using previously acquired transcriptomics datasets (Supplementary Table 4). Transcriptomics datasets were processed using the same methods described for Xenopus transcriptomics and normalized using a series of transformations. Because RNA-seq data typically have a wider dynamic range than microarray datasets, RNA-seq and microarray data was transformed into a common space²⁶. RNA-seq datasets were transformed using a hyperbolic sine function and the log₂ transformation was applied to microarray data²⁶. Subsequently, a variance stabilizing transformation was applied to remove potential mean-variance relationships in each dataset²⁷. Common genes expressed and measured between datasets (926 genes) were considered for further analysis. Dimensional reduction was performed using Principle Component Analysis (PCA) and visualized using Uniform Manifold Approximation and Projection (UMAP) in Python with the help of sklearn and UMAP packages.

Commonalities in gene expression values across species were analyzed by plotting heatmaps of differential gene expression using the Python packages seaborn and matplotlib. All datasets were normalized as described above and log fold change for each gene was calculated by combining samples of a condition using geometric mean. The Mann-Whitney U-test was used to compute the p values of differential expression, which were FDR-corrected. To construct a common heatmap across datasets, rows missing data (NaN,Not-a-Number) were dropped from all differential gene expression datasets. The log fold change column from differential genes of all species also were clubbed together, and the variance in log fold change was calculated for all genes. The genes (96 genes) with variance greater than the mean variance (in fold changes) were used to construct the heatmaps.

Transcriptomics gene networks were constructed using GENIE3. All cell lines were used to investigate potential changes across the entire gene network, regulator genes were not down selected. NeMoCAD takes on an agnostic approach and does not down select any cell types. Adjacency matrices were visualized for each network and the distribution of edge weights in the network was assessed by histogram and boxplots in MATLAB. Digraph objects were constructed from adjacency matrices and visualized in MATLAB for the strongest connections in each network (e_w > 99.99% for each network). Resultant networks were combined across all datasets for each species and plotted with different colors designating each dataset. A cross-species network that includes the strongest network connections across all datasets was formed by combining each species-specific dataset and plotted with different colors designating each species. To identify the strongest regulator and target gene nodes in the resultant cross-species network, the in-degree and out-degree were calculated in MATLAB. The in-degree is defined as the number of edges with that node as the target and out-degree of a node is equal to the number of edges with that node as the source.

Drug screening

Drugs used for tadpole screening including vorinostat, trofinetide, ivermectin, and clozapine were purchased from Sigma-Aldrich (St. Louis, MO, USA) or MedChem Express (Monmouth Junction, NJ, USA) and dissolved in DMSO to a stock concentration of 100–1000 mM, depending on the solubility. We leveraged our previous experience performing high-throughput tadpole screenings²⁸ to determine the concentration range to test for efficacy and toxicity. Compounds were diluted to their final concentrations with a 0.1% DMSO concentration. Control groups included MeCP2 KD and wild-type age-matched tadpoles that were treated with the test compound and vehicle (0.1% DMSO concentration). Healthy and infection free tadpoles were randomly selected for each dosing group. Media with dosed drug were made fresh and exchanged every 2 days following feeding. Dosing began at stage 45 and lasted for 7 days, after which standard 0.1X MMR was used to evaluate drug washout effects. Free-swimming tadpoles were imaged in 60 mm dishes using a SONY Alpha a6100 camera with 16 mm objective against an illuminated background. Tadpole behavior was scored manually after a 10 min acclimatization to the illuminated background at 15 °C by noting, within a 20 min period, the most severe seizure-related phenotype based on the PTZ induced seizure tadpole model that develop class C- shaped contraction¹⁴. Researchers were blinded to the condition in each video and scored based on key morphological signatures.

X. laevis tissue processing

After stage 46 tadpoles were euthanized with 20x Tricaine (3.36 g/L), tadpoles were washed three times with PBS (-/-) and fixed in 4% Paraformaldehyde in PBS + Mg + EGTA (MEMFA) (Alfa Aesar, Tewksbury, MA, USA) for 2 h at RT. For sections, Xenopus tadpoles were cryoprotected by sequential incubation of 10%, 20% and 30% sucrose in 1× PBS(-/-), only transferring to the subsequent concentration when the samples sunk to the bottom of the scintillation vial. Cryoprotected tadpoles were embedded in 5% agarose and were stored in −80 °C. For GI staining, the gastrointestinal tract was microdissected with a stereoscope before embedding the tissue, which were subsequently sectioned coronally at 20 μm thickness on a Cryostat (Leica, CM3050 S, Wetzlar, Germany) and mounted on superfrost plus slides (Thermo Fisher, Waltham, MA, USA). Slides were dehydrated for 30 min at 20 °C, before storing at −80 °C.

In situ hybridization

Following X. laevis tissue processing of cryosections, sample pretreatment, RNAscope Target Retrieval and the RNAscope Assay was conducted by the HMS Neurobiology Imaging Facility following manufacturer’s instructions. Xl.MeCP2 was custom designed by Advanced Cell Diagnostics. Rpl8 was used as a positive control (Xl-LOC108706872-C3, Cat No. 516501-C3, Advanced Cell Diagnostics, Hayward, CA) and DapB as a negative control (Cat No. 320758, Advanced Cell Diagnostics, Hayward, CA).

Immunohistochemistry of Xenopus laevis sections and whole mounts

We utilized immunohistochemical methods to characterize the Rett tadpole model postmortem as live confocal imaging was not possible due to the high frequency of seizures. Following tissue processing, cryosections were blocked for 2 h at room temperature and subsequently incubated with primary antibodies overnight at 4 °C, including mouse anti-alpha acetylated tubulin (1:1000 dilution, Sigma Aldrich, Burlington, MA USA), mouse anti-beta-tubulin (1:1000 dilution, E7, Developmental Studies Hybridoma Bank, Iowa City, Iowa, USA), anti-MeCP2 (1:1000 dilution, Invitrogen, Waltham, MA USA), Hoechst (1:1000 dilution, ThermoFisher, Waltham, MA, USA), negative control anti-mouse IgG (1:1000 dilution, Invitrogen, Waltham, MA, USA), and anti-rabbit (1:1000 dilution, Invitrogen, Waltham, MA USA), as well as FITC-conjugated isolectin B₄ (1:1000 dilution, Thermo Fisher, Waltham, MA, USA). Slides were washed twice with TRIS-buffered saline, with 1% Tween-20 (TBST) for 5 min each and incubated in the dark with their corresponding secondary antibody: Goat anti-rabbit Alexa 647 (1:1000 dilution, Abcam, Waltham, MA, USA) and Goat anti-mouse Alexa 594 (1:1000 dilution, Abcam, Waltham, MA, USA) for 1 h at room temperature. Slides were washed three times with TBST for 5 min each on an orbital shaker and mounted in ProLong Gold Antifade Mountant (ThermoFisher, Waltham, MA, USA) with a #1 coverslip and sealed. Slides were stored in the dark at 4 °C before image acquisition. Z-stacks were acquired at 1 μm to visualize brain and gastrointestinal sections.

X. laevis embryos were processed for whole mount immunofluorescence after fixing the tadpoles in 4% paraformaldehyde in PBS + Mg + EGTA (MEMFA) (Alfa Aesar, Tewksbury, MA, USA) for 1 h at RT. Subsequently, fixed tadpoles were quenched for 1 h in a lightbox in 5% formamide and 4% hydrogen peroxide in PBS, permeabilized for 1 h in PBS(-/-) with 0.1% Triton X-100 (PBST) and blocked for 1 h in 10% CAS-Block in PBST (Life Technologies, Carlsbad, CA, USA)²⁹. Primary antibodies were incubated overnight at 4 °C, including mouse anti-α acetylated tubulin (1:700 dilution Sigma Aldrich, Burlington, MA USA), mouse anti-β-tubulin (1:100 dilution, E7, Developmental Studies Hybridoma Bank, Iowa City, Iowa, USA), Hoechst (1:100 dilution, ThermoFisher, Waltham, MA, USA) and anti-mouse IgG (1:100 dilution). Slides were washed in PBST, blocked in 10% CAS-Block for 30 min, and incubated with the corresponding secondary antibody for 2 h at room temperature including Goat anti-mouse Alexa 546 (1:100 dilution, ThermoFisher, Waltham, MA, USA) and Goat anti-rabbit 633 (1:100 dilution, ThermoFisher, Waltham, MA, USA). Tadpoles are washed with PBST for 1 h and then washed in PBS for another hour before mounting in Vectashield (Vector Laboratories, Burlingame, CA, USA). Z-stacks were acquired at 1 μm for visualizing epidermal MCC’s and 3 μm for visualizing entire brain regions.

Drug target discovery in X. laevis using thermal proteome profiling

Sample preparation

While we largely followed the protein integral solubility alteration (PISA) thermal proteome profiling protocol³⁰, we adapted the mammalian cell- and tissue-specific parameters to tadpoles by optimizing the method to tadpoles by reducing thermal treatment temperatures to account for the lower normal temperature range of X. laevis. Unmodified X. laevis tadpoles (stage 47–50) were exposed to 25 mM vorinostat or 0.1 mM ivermectin for 2 h while freely swimming in 0.1x MMR medium. Following euthanasia using 20x Tricaine solution (containing 3.36 g/L of Tricaine buffered in 1 M Tris to pH 7), animals were washed with fresh 1x MMR, cut into thirds (head, abdomen containing all visceral organs, and tail), and placed into individual Eppendorf tubes covered with MMR containing the same drug and concentration and incubated for 3 min at one of 9 temperatures spanning 30–60 °C. 30 °C was selected as the control temperature, where no appreciable thermal denaturation should take place, based on the maximum tolerable temperature for X. laevis³¹. Following aspiration of the remaining liquid, samples were flash frozen in liquid nitrogen.

Samples were further processed and analyzed by Phenoswitch Bioscence, Inc. (Sherbrooke, Québec, Canada). Thermally-treated samples were lysed in 100 µL PBS with 0.4% NP-40 with 4 cycles of freeze/thaw. Insoluble material was cleared by centrifugation (10 min, 10,000 G, 4 °C). 11 µL of supernatant for each of the 9 temperatures were then pooled in one tube and samples were centrifuged again at 13,000 RPM for 75 min, at 4 °C to pellet precipitated proteins. 80 µL of the supernatant were reduced with 10 mM DTT for 15 min at 65 °C and alkylated with 15 mM IAA and 30 min at room temperature in the dark. Proteins were precipitated with 8 volumes of ice-cold acetone and 1 volume of ice-cold methanol overnight. Protein pellets were washed 3 times with 250 µl of ice-cold methanol and resuspended in 100 µL digestion buffer. Digestion was carried for 4 h in 50 mM Tris ph 8 + 0.75 mM Urea + 1 µg trypsin/LysC at 37 °C with agitation. Another 1 µg of trypsin/LysC was added, and digestion was continued overnight. Peptides were purified by reversed phase SPE and analyzed by LC-MS/MS.

Mass spectrometry

Acquisition was performed with a TripleTOF 6600 (Sciex, Foster City, CA, USA) equipped with an electrospray interface with a 25 μm iD capillary and coupled to a Micro LC200 (Eksigent, Redwood City, CA, USA). Analyst TF 1.8 software was used to control the instrument. Acquisition was performed in data independent acquisition (DIA or SWATH) using gas phase fractionation (GPF1 from 350 m/z to 800 m/z and GPF2 from 800 m/z to 1250 m/z). The source voltage was set to 5.5 kV and maintained at 325 °C, curtain gas was set at 45 psi, gas one at 25 psi and gas two at 25 psi. Separation was performed on a reversed phase Kinetex XB column 0.3 μm i.d., 2.6 μm particles, 150 mm long (Phenomenex) which was maintained at 60 °C. Samples were injected by loop overfilling into a 5 μL loop. For the 60 min LC gradient, the mobile phase consisted of the following solvent A (0.2% v/v formic acid and 3% DMSO v/v in water) and solvent B (0.2% v/v formic acid and 3% DMSO in EtOH) at a flow rate of 3 μL/min. Both GPF files were analyzed on a previously-generated 3D ion library using the SWATH 2.0 microapp from Peakview (Sciex, Foster City, CA, USA). Each GPF file was analyzed using 10 peptides per protein, 4 MS/MS transition per peptide, 12.5 min RT window and 25 ppm XIC width. The reported quantification for a protein is the sum of all the correctly integrated peptides (FDR < 0.05) in both GPF files.

Mouse studies

Mouse care

All animal experiments and procedures were reviewed and received ethical approval by the HMS Institutional Animal Care and Use Committee regulations (Protocol # IS00002714). Every effort was made to minimize their suffering. Experiments were performed with exon 3–4 knockdown MeCP2-null (MeCP2^-/Y) (JAX strain 003890) and age-matched WT male littermates (MeCP2^+/Y), purchased from Jackson Laboratories (Bar Harbor, ME). We chose this commonly used preclinical model due to its well-defined disease onset and phenotype. Since MeCP2-null female mice are rare (<1% of all offspring)³², we performed this study only on male mice. Future studies could include using female hemizygous MeCP2^-/X mice to address other aspects of the disease (e.g., sex and mosaicism). All mice were housed in ventilated racks under specific-pathogen-free conditions at a room temperature in a 12/12 h light/dark cycle with food and water ad libitum. Mice were assessed for disease severity and endpoint criteria on a daily basis by weighing and using the phenotypic scoring method developed by the Bird laboratory^10,33, which includes assessing mobility, gait, breathing, tremor, and general condition. The only exclusion was the humane endpoint, which was determined by 20% body mass loss or a score of 2 in criteria D, E, or F. Scores were produced by 1 blinded researcher and 1 non-blinded researcher. Because of this and the fact that the MeCP2 KD were very different from wild-type controls, care was taken during data analysis to use objective metrics whenever possible to avoid bias. All mouse studies were repeated two or more times, except the mouse study with oral dosing, which was conducted once. Animals were acclimatized for 2 days before beginning dosing. Dosing and Behavioral procedures began when mice reached 4 weeks of age, and where cages of mice were selected at random for each condition using = rand() with excel. Individual mice were identified over time by marking their tails with a permanent marker.

Drug efficacy testing

Vorinostat and trofinetide for mouse studies were purchased from MedChem Express (Monmouth Junction, NJ, USA). We utilized mouse to human dose conversions and doses previously used in studies of vorinostat in mouse models to determine the dose and concentration range to test for therapeutic efficacy. For intraperitoneal injection, trofinetide was solubilized in DMSO, and vorinostat was solubilized in a method using 2-hydroxypropyl-β-cyclodextrin powder (HPβCD) was purchased from Acros Organics (ThermoFisher, Waltham, MA, USA)^34,35. Vorinostat was dissolved in 100 mM HPβCD by boiling the solution at 95 °C for 5 min. Mice were dosed intraperitoneally at 100 mg/kg with trofinetide and at 100 mg/kg and 50 mg/kg with vorinostat on a daily basis. For oral dosing, commercial instructions were followed to add solubilized vorinostat in HPβCD (100 mg/kg) to Medigel (purchased from ClearH20, Westbrook, ME, USA) and were replaced every other day. Wild-type littermate gender matched controls groups were included to test the vehicle and the compounds (Vorinostat and Trofinetide). A group of MeCP2 KD male mice were also treated with vehicle controls (0.1% DMSO) and 100 mM HPβCD for trofinetide and vorinostat, respectively.

Behavioral testing

We performed a battery of behavioral tests at the Harvard Medical School Behavioral Core facility to capture the behavioral paradigms observed in Rett patients. There were a few tests that were not compatible with the MeCP2-null’s smaller stature, gait modifications and our dosing regimen. For example, to assess sociability, we attempted the 3-Chamber Test, which assesses a normal mouse’s tendency to spend time with a non-littermate mouse, the stranger mouse, contained in a metal wire cup. This is performed after habituating the test mouse with the metal cup alone in the test chamber. We found it difficult to determine a clear phenotype with the MeCP2 null mice because of its tendency to climb above the metal cup, which allowed it to avoid the stranger mouse despite its close proximity. This is something that we did not observe in age and gender matched wild-type controls (Supplementary Fig. 8f). Thus, the tests shared below produced the most consistent and robust phenotype to enable us to determine therapeutic efficacy. For each behavioral test, the entire mouse cohort was tested within a 4-h period and individual mice were randomized in the order of testing. During each experimental run, there were instances where data could not be collected by the entire group either when mice could not be identified by the software for tracking or when individual mice failed the positive controls or assumptions of the test. To avoid habitation bias, a failed test could not be repeated.

Elevated plus maze (EPM)

The elevated plus maze (EPM) consisted of two open and two closed (34 cm long and 5 cm wide) arms extended out from a central platform 50 cm above the floor. The test was carried out in dim ambient lighting. The room of the experiment and its spatial cues were kept consistent and minimal. Mice were habituated in the test room for 30 min before the start of the test. Mice were placed near the center compartment of the maze, facing an open arm, and allowed to explore the apparatus for 5 min. After each test, the maze was cleaned thoroughly with 70% ethanol and left to dry before the start of the next test to remove odorant cues. A computer-assisted video-tracking system (EthoVision XT 14, Noldus, Leesburg, VA) was used to record the number of open and closed arm entries as well as the total time spent in open, closed, and center compartments. An increase in the percent time spent or entries into the open arms was used as a surrogate measure of anxiolytic-like behavior³⁶. Data could not be collected from mice that fell off the arms.

Spatial novelty y-maze (Y-maze)

The spatial novelty y-maze consisted of a Y-shaped maze with three arms at a 120° from each other as well as a removable blockade used to prevent access to one of the arms during the Habitutation phase. The test was carried out in dim ambient lighting. The room of the experiment and its spatial cues were kept consistent. One of the three arms was defined as the start arm, where the mice would be placed at the start of the experiment and was kept consistent throughout the entire study. The test consisted of a 3-min Habituation phase where one of the two non-start arms was blocked off, a 2-min intertrial interval (ITI) outside of the maze and a 3-min Test phase, where the previously blocked arm was exposed to the mouse. The arm that was blocked off during the Habituation phase was randomized and noted for each mouse.

Mice were kept in an adjacent room to not disturb the test. The test mouse was placed in a holding cage and brought to the test room with the maze and immediately placed into the Y-maze at the designated start arm and allowed to explore the start and familiar arm trial for 3 min. At the end of the habituation phase, the mouse was placed into the holding cage for 2 min ITI. During this time, the blockade was removed, and the maze was cleaned with an ammonia-based cleaner with a paper towel and left to dry to remove odorant cues. The mice were then placed back in the start arm and left to explore the entire maze for 3 min. After the test trial, the mouse was transported back to its home cages and the maze was again wiped dry with an ammonia based cleaner. A computer-assisted video tracking system (EthoVision XT 14, Noldus, Leesburg, VA) was used to record the number of entries, time spent in each arm and as distance traveled in each arm during the Habituation and Test Phase. The percentage of time and distance spent in the novel or previously blocked arm was used to assess spatial novelty seeking behavior³⁷.

Mouse histology and immunostaining

Mice were euthanized with CO₂ inhalation, transcardially perfused with PBS (-/-) and 4% PFA, and tissues were removed and fixed with 4% paraformaldehyde overnight at 4 °C and stored in 100% ethanol at 4 °C before paraffin processing. Regions of the gastrointestinal tract and lung were embedded into paraffin blocks by the Beth Israel Deaconess Medical Center’s Histopathology Core Facility^38,39,40. Paraffin blocks were sectioned into 15 μm thickness. Mouse lung tissues were sectioned coronally and GI tissues were stained with hematoxylin and eosin (H&E) by the Beth Israel Deaconess Medical Center’s IHC Core Facility. For immunofluorescence, the slides were stained for rabbit anti-acetylated α-tubulin 1:1000 (Abclonal, ab179484, Woburn, MA, USA) mouse anti-CD-64 1:1000 (R&D systems, AF3628, Minneapolis, MN, USA), chicken anti-βIII-tubulin 1:1000 (Biolegend, 801202, San Diego, CA, USA) immunofluorescence by Beth Israel Deaconess Medical Center’s IHC Core Facility.

Mouse brains were fixed with 4% paraformaldehyde at 4 °C overnight, cryoprotected in a sucrose gradient (10%, 20% and 30%), and subsequently embedded in OCT Tissue Tek. The embedded tissues were subsequently sectioned coronally at 30 μm and the free-floating sections were blocked with 10% CAS Block in PBST for 1 h, stained with primary antibody (rabbit anti-iba-1, 1:1000, Wako, #019-19741; mouse anti-NeuN, 1:1000, Abcam, #ab104224; Hoechst, 1:1000, ThermoFisher, H1399) overnight at 4 °C. After three 10-min washes in PBS (-/-) on an orbital shaker, the sections were stained in secondary antibody (anti-rabbit Alexa 647 for iba1, anti-mouse Alexa 594 for NeuN, at 1:1000) for 1 h at room temperature. After three 10-min washes on an orbital shaker, the sections were mounted on superfrost plus slides (ThermoFisher, Waltham, MA, USA) in ProLong Gold Antifade Mountant (ThermoFisher, Waltham, MA, USA) with a #1 coverslip and sealed. Z-stack images were acquired at 1 μm intervals using the 63× objective with glycerol⁴¹.

Multiplex chemokine assay

After mouse euthanization, blood was collected by cardiac puncture using at 23G needle in ethylendiamine tetra-acetic acid (EDTA) treated tubes (Microvette 200 K3EDTA, CAT# 20.1288.100, Sarstedt, Numbrecht, Germany). To avoid hemolysis, the tubes were kept at room temperature and quickly processed for plasma. Tubes were centrifuged at 2000 × g for 10 min at RT and the plasma was separated and stored in −80C until assayed. Chemokine levels were measured in mouse plasma using the Bio-Plex Pro Mouse Chemokine Panel 31-Plex (Cat# 12009159, Bio Rad, Hercules, CA, USA). Measurements were performed using the Bio-Plex 3D Suspension Array System (BioRad, Hercules, CA, USA), following the manufacturer’s instructions.

Microscopy and image processing

H&E slides were imaged with 4× and 20× objectives with the Biotek Cytation 5. Immunofluorescent Xenopus whole mounts and Xenopus and mouse-stained sections were imaged with the Leica SP5 X MP Inverted Confocal Microscope using 25× (NA 0.95) water objective and 63x (NA 1.3) Glycerol objectives. Sections were imaged in white light and diode laser scanning mode and whole mounts were imaged using multiphoton pulsed IR laser scanning. Serial scanning was used during acquisition to avoid bleed-through. Image settings were kept constant between groups, in terms of laser power, bit-depth, z-stack interval, resolution (pixels), aperture settings, and gain/offset. Images were post-processed and analyzed using ImageJ (Fiji, NIH).

MCC morphological analysis

Multiple 2-photon z-stack X. laevis epidermis MCC images were acquired at 25× and 63× and identified with acetylated α-tubulin staining, where the apical surface and the cilia were identified. Z-stacks where vertical cross sections of each MCC were identified, collapsed at maximum intensity, background subtracted, and subsequently assessed for cilia orientation and cilia length relative to the basal body⁴². After using a minimum binding box algorithm, the polygon vertices were utilized to trace the basal body and cilia for each MCC. These traces were subsequently used to measure the length of each cilium defined by the apical vertices to the basal body and its orientation by measuring the angle formed from the cilium relative to the basal body axes for each MCC.

MCC functional assay

We visualized the cilia-driven flow of fluorescent microspheres over the surface of the xenopus tadpole by taking times series with the Leica SP5 x MP inverted confocal. Our methods were based on previous reports⁴³ with some modifications. Stage 45 tadpoles were placed in a 24 well μ-Plate 14 mm (ibidi, Gräfelfing, Germany), and underwent anesthesia by incubation with 1× tricaine (0.168 g/L) before proceeding to image. One tadpole was placed per well for individual tracking. Fluorescent microspheres beads (Thermo Fisher, Waltham, MA) were added at a final dilution of 1:100, at 0.02% solids. We chose the bead size of 3 μm, because during preliminary testing and as is reflected in the data, WT X. laevis tadpoles consistently was able to clear beads that were 3 μm diameter at this concentration. The Leica SP5 × MP laser scanning mode was used to visualize the fluorescent microspheres beads at an excitation of 520 nm and emission of 560 nm with the following settings: 512 × 512 pixels, bidirectional scanning mode, and at a scan speed of 800–1600 Hz, with the 25× (NA 0.95) water objective, and over the time course of 5 min every 5 s. Time series were acquired of in the olfactory and lateral region of the epidermis.

Quantification of bead displacement

The time series stacks were collapsed at maximum intensity to visualize the overall bead flow. Intensity was measured as a function from the distance from the surface of the epidermis.

Mean intensity quantitation

To generate a ROI mask, channels acquired in one z-stack were split and the maximum intensity projection of a ‘control’ channel was taken. In the case for assessing αaT expression mouse GI tract and bronchiole, the βIII-tubulin channel was used, and for assessing MeCP2 expression, the \(\alpha\) -acetylated tubulin channel was used in the image stack. Subsequently, background noise was subtracted, and auto-thresholded, whose settings were held constant between groups. The ROIs generated from the control channel were used to quantify the maximum projected channel of interest, to obtain the mean intensity (integrated density per area). A minimum of 3 serial sections were analyzed per animal.

Counting Ib4+ cells

Confocal Z-stacks of X. laevis GI tract images were quantified for Isolectin B4 (ib4)+ cells, with the maximum intensity projection of the ib4+ channel. Following noise reduction, and auto-thresholding, ib4+ cells were segmented and counted using the particle analysis function in the ImageJ software (NIH, Bethesda, MD).

Sholl analysis

Confocal (at 1 μm intervals) Z-stacks were collapsed at the maximum intensity in the iba-1 channel, noise de-speckled, subsequently segmented using the auto-threshold function, and microglia were duplicated into 8-bit TIF files, in the field of view, ensuring that for each mouse at least 10 microglia were selected in the ipsilateral and contralateral region of the olfactory bulb. Using the Sholl analysis plugin in ImageJ (NIH, Bethesda, MD), traces and intersections were generated.

Statistical analysis and reproducibility

Statistical analysis of biological results was performed using Graphpad Prism 9.3. For comparison of two conditions, Student’s t test was used unless otherwise indicated and F-test was used to compare variances. ANOVA with pairwise t-tests and Bonferroni correction or the Holm-Šídák tests was used to evaluate significance of multiple variables. P values < 0.05 were considered significant. Tadpole screens and mouse studies were conducted with a minimum of 5 animals each and repeated except for the mouse oral dosing study to account for animals that might reach humane endpoints or if data could not be collected during the behavioral study,. All data indicate separate biological replicates (i.e., tadpole or mouse) unless otherwise noted, where a minimum of N = 3 animals were used, with N ≥ 3 technical replicates each (e.g., histology sections).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Modeling Rett syndrome genetic and phenotypic heterogeneity in Xenopus laevis

A major challenge to Rett therapeutics discovery lies in the inherent heterogeneity observed in patients stemming from gene expression variability due to the gene regulatory impact of distinct MeCP2 mutations, combined with tissue heterogeneity caused by mosaicism resulting from X chromosome inactivation. As a result, Rett patients exhibit a spectrum of disabilities. To establish an animal model of Rett that would incorporate the heterogeneity observed in patients, we used CRISPR to generate a mosaic knockdown of MeCP2 protein in Xenopus laevis tadpoles, which have been previously used to model neurodevelopmental diseases due to their evolutionary proximity to mammals and well-characterized neural development^44,45.

Injection of X. laevis embryos at the 4-cell stages using cas9 pre-complexed with CRISPR sgRNA targeting 6 sites of the MeCP2 gene, including the regions corresponding to the DNA-binding ___domain, resulted in a tunable model of MeCP2 knockdown with a heterogeneous phenotype when assessed at the population level. Based on our target selection and injection strategy at the 4-cell stage, each individual animal was edited as a mosaic and presented varying degree of severity, without apparent toxicity from Cas9 or sgRNA. Indeed, all embryos developed to the swimming tadpole stages (Nieuwkoop-Faber stages 45–50) with no apparent loss of viability or morphological defect in earlier stages.

Importantly, swimming tadpoles with MeCP2 knock down (Rett tadpoles) exhibited a broad range of abnormal behavior compared to wild type vehicle-injected control tadpoles (Control) including darting motions and rapid repetitive swimming in tight circles (Fig. 1a, b and Supplementary Movie 1) reminiscent of the repetitive motions observed in human Rett patients as well as C-shaped or straight rigor (Supplementary Movie 1) that has been shown to correspond to seizures in prior Xenopus studies⁴⁶. The cas9 protein and sgRNA doses were titered to achieve between 5 and 95% gene editing efficiency as determined by qPCR without resulting in toxicity (Supplementary Fig. 2b). All tadpoles exhibited MeCP2 gene insertion/deletion polymorphisms (indels) at each of the MeCP2 gene target sites (Supplementary Fig. 2c) and reduced MeCP2 RNA expression (Fig. 1c). We utilized immunohistochemical staining to quantify levels of MeCP2 protein levels across tissues and to address the mosaicism resulting from the stochastic gene editing, where MeCP2 protein levels were also heterogeneously and significantly reduced in the brain (Fig. 1d, e) and olfactory multiciliated cells (MCC) (Fig. 1d, f). Thus, this editing approach resulted in an animal model of Rett with intra-subject and inter-subject diversity at molecular and phenotypic scales similar to that observed in humans.

a, b Tadpole models of Rett syndrome exhibit distinct swimming behavior in 60 mm diameter dishes following MeCP2 knockdown (Rett) compared to buffer-injected controls (Control). MeCP2 RNA expression in brain and gastrointestinal (GI) tract using RNAscope (tadpoles per condition: N = 5 Rett, N = 6 Controls) c and MeCP2 protein in brain (tadpoles per condition in the forebrain: N = 5 Rett, N = 5 Controls; in midbrain: N = 5 Rett, N = 8 Controls) d, e and multi-ciliated olfactory cells (tadpoles per condition: N = 5 Rett, N = 4 Controls) f using immunohistochemistry (rabbit anti-MeCP2 and goat anti-rabbit Alexa Fluor 647), are significantly reduced in Rett syndrome tadpole models while maintaining a large degree of heterogeneity spatially d and across tissues, e, f; scale bars, 100 µm in brain, 20 µm in olfactory multi-ciliated cells. All error bars indicate standard deviation g Volcano plot of differentially expressed genes showing that only 70 genes are significantly up- or down-regulated following MeCP2 knockdown (p < 0.05, fold change > 2, t-test with Bonferroni correction, N = 3). Gene regulatory networks for control, h, and Rett i, tadpoles reveal large network rearrangements with an increase in BDNF centrality. j Comparison against other Rett syndrome models and clinical samples indicates minimal differential gene expression across the 96-gene signature developed to classify Rett syndrome. k Network-level comparison across species and tissues (left) and identified shared target (middle) and regulator (right) genes.

Full size image

In line with prior studies in humans and other model organisms⁴⁷, when we analyzed the transcriptome-wide effects of MeCP2 knockdown on the Rett phenotype in Xenopus, only subtle shifts in gene expression were detected with 70 out of 10,935 genes probed undergoing a significant change in gene expression (p_adjusted< 0.05) (Fig. 1g, Supplementary Table 3). Of the 70 differentially expressed genes, 37 are involved in metabolic processes, 9 operate within developmental processes, and 8 are part of signal transduction pathways (Supplementary Fig. 3a). This aligns with changes in genes that control metabolism or regulate neuronal processes (e.g., ion transport, nervous system development) in Rett patients².

Interestingly, analysis of gene co-expression networks identified substantial reorganization from control animals (Fig. 1h) following MeCP2 knock down (Fig. 1i), characterized by a loss of strong connectivity among genes involved in the regulation of MeCP2 (FMR1 and TRKB receptor NTRK2^48,49), neuronal development (BDNF ^50,51), and circadian rhythm (PER1, PER2⁵²) all of which are altered in Rett patients^50,53. Notably, after MeCP2 knock down, BDNF strengthened as a network hub and exhibited greater co-expression connectivity amongst gene nodes (Fig. 1i).This was unexpected because although BDNF is involved in neuronal development, synaptic transmission, and plasticity, it declines with the onset of Rett-like neuropathological and behavioral phenotypes⁵⁰. However, BDNF is a known target of repression by MeCP2 and our findings are consistent with prior work suggesting that downregulation of BDNF is a later and indirect outcome of MeCP2 deficiency^50,51.

We then comprehensively compared the effects of MeCP2 knock down in Xenopus to changes observed in published transcriptomic data sets from Rett patients as well as MeCP2 knockout mouse and zebrafish models (Supplementary Table 5). After performing dimensionality reduction in this cross-species analysis, K-Means clustering applied to the principal components revealed 3 clusters, with mouse and human clustering together and zebrafish and Xenopus forming 2 independent clusters (Supplementary Fig. 3b). Importantly, diverse gene-expression abnormalities that occur across cell types in Rett syndrome could contribute to the differences observed across these datasets². In this case, both the human and mouse transcriptomics were performed on brain tissues, whereas the zebrafish and Xenopus data are from whole organisms, which may account for the clustering patterns observed.

Despite evidence of species-specific clustering, a common gene signature comprised of 96 genes was identified across all species by searching for genes with a fold change variance greater than mean variance. However, the common genes within this transcriptomic signature are expressed at similar low fold change levels in 7 of the 9 datasets analyzed, with slightly greater variations in expression only observed in 1 out of 5 mouse datasets and 1 out of 3 human datasets (Fig. 1j). This finding provides additional evidence that mutation of MeCP2 produces widespread, yet subtle, perturbations in gene expression.

Using gene co-expression network analysis²⁵, we also identified interconnected gene sub-networks that are conserved across species and tissues in the MeCP2 deficient condition (Fig. 1k, Supplementary Tables 4 and 5). Comparison across all species identified several shared major gene targets and regulators despite differences in species-specific network structure (Supplementary Fig. 3c). Interestingly, WNT8B, which is known to be modulated by the Rett-implicated gene Foxg1⁵⁴, was found to be the most interconnected gene target within the MeCP2 cross-species network. The strongest gene regulators in this network, NR4A3 and SOD2, are both involved in oxidative metabolism, which is also known to be dysregulated in Rett Syndrome⁵⁵. Thus, this CRISPR-enabled Xenopus model of Rett replicates many genetic and phenotypic features of the disease seen in humans.

Computational discovery of drugs that reverse the network-level impact of Rett

In Rett, MeCP2 knockout results in widespread but subtle gene expression patterns, and most effects are not detectable using standard statistical analysis. Researchers have applied network-based analytical methods to uncover previously-hidden connections within RNA and protein interaction networks^47,53,56, but these approaches have not been applied to Rett syndrome therapeutics discovery. To leverage transcriptomics analysis for drug discovery, used a computational framework we termed nemoCAD (network model for causality-aware discovery)(Fig. 2a) that compares transcriptomic signatures of disease versus healthy control subjects¹⁵ to predict compounds from the Library of Integrated Network-Based Cellular Signatures (LINCS) database⁵⁷, including FDA-approved drugs, which have a high likelihood to reverse the disease signature state back to a healthy state. Though signature-based tools have been previously developed to enable drug discovery (e.g., CMAP⁵⁸, NicheNet⁵⁹, GPSnet⁶⁰), their major drawback is the impact of the expression of any single gene or protein on the resulting prediction. This is because current computational approaches rely on known target genes and gene signatures of a perturbation without abstraction, which can make such tools susceptible to bias and limited in drug coverage^61,62. While the adoption of Weighted Gene Coexpression Analysis aims to improve such existing tools for intramodular connectivity, it continues to rely upon the presence or absence of featured target genes, and hence is still restricted to an undirected network⁶³. Network-based methodologies have helped to address these limitations of relying on the sole presence or absence of a gene⁶⁴ which underscores the importance of a network-based abstraction. Furthermore, to be used effectively, these tools require well-matched disease datasets to the human cell line-based datasets present in LINCS, which may not be feasible for many diseases.

a Network model for causality-aware discovery (nemoCAD) combines a directed gene-gene and drug-gene interaction network extracted from CTD, TRRUST, and KEGG databases with interaction probabilities inferred from single gene and drug perturbations in LINCS. Transcriptome data from any disease model or patient and corresponding control are used to identify the relevant subnetwork and disease-specific node weights that account for probabilities of up-/down-regulation of a gene. A drug-gene interaction probability matrix, inferred from LINCS, is computationally screened against the disease-specific subnetwork to identify compounds that significantly interact with the subnetwork and are ranked by their predicted ability to restore the disease transcriptome back to a healthy state based on single gene and gene network signatures. Downstream analyses can be performed on the resulting gene-gene interaction subnetwork by interrogating the underlying subnetwork structure to find control nodes and other network metrics. Additionally, the chemical structures of the predicted drugs can be clustered by structural similarity based on SMILES notation and annotated protein targets and pathways from DrugBank data. Created in BioRender. Lin, T. (2025) https://BioRender.com/40iq3jv. b Graph showing relative effects on seizure score over 10 days of treatment in MeCP2 KD and buffer-injected tadpoles of vehicle and 25 µM vorinostat (Rett-Vor) versus 70 mg/mL trofinetide (Rett-Tro), a clinical-stage drug with demonstrated efficacy (tadpoles per condition: N = 13 Rett-Vor, N = 20 Rett-Vehicle, N = 8 Rett-Tro). Vehicle-treated MeCP2 KD tadpoles did not survive past day 3 of the treatment period in one study and past day 6 in a second study. Wild-type Xenopus buffer-injected Controls treated with vehicle and vorinostat did not exhibit seizures. c Changes in swimming score at day 6 of treatment vs. baseline (day 0) is significantly improved by both drugs in Rett tadpoles. Error bars indicate standard deviation.

Full size image

In contrast, nemoCAD pursues an entirely different approach to transcriptomic signature analysis that is based on a Bayesian network encoding directed gene-gene and drug-gene interactions combined with machine learning¹⁵. This approach simulates network impacts of each drug, which covers any gaps in the network. By leveraging data in LINCS containing transcriptomics results following single gene knockout/knockdown and overexpression as well as addition of drug perturbagens in human cells, we calculated the probabilities of pairwise interactions among all of the genes and drugs in a human gene-gene and drug-gene interaction network whose structure was extracted from Comparative Toxicogenomics Database (CTD)⁶⁵, KEGG⁶⁶, and TRRUST⁶⁷. Unlike signature based computational tools that take gene expression signatures at face value, this probabilistic approach offers greater abstraction while prioritizing for the network impact of perturbations, including those caused by disease or addition of drugs, and thereby reduces the impact of tissue-specific or species-specific transcriptomic biases as well as the influence of the expression of any individual gene. The NeMoCAD network approach is resistant to noise, which can stem from technical or biological sources, as we saw no deterioration in drug prediction performance up to 50% added stochastic noise (Supplementary Fig. S1).

To identify drugs that might reverse disease states, we generated probabilistic network maps for each drug in the aggregate 19,800+ compound dataset. By inputting and comparing gene microarray data from MeCP2 knockdown tadpoles versus vehicle-injected controls, we identified sets of differentially expressed genes contained within the gene-gene interaction network with a p < 0.05 and a log2 fold change between 0.7 and 3.1 in 0.3-fold-change increments. For each fold-change, we generated a ranked list of compounds predicted to induce reversal of transcriptome-wide changes of: (1) single gene states using a cross-correlation score, (2) single gene states using a cross-entropy score, and (3) an overall gene network state score corresponding to the probability of involvement of the drug nodes in the identified subnetwork. These drug lists were combined to identify drugs that were ranked by being most consistently predicted to reverse the Rett transcriptome-wide changes and robust across the computational settings (Supplementary Table 6). The prediction of vorinostat was further benchmarked against predictions from CMAP (Supplementary Fig. S1) and the results indicated that this drug would not have been highly predicted using that approach. Drugs were manually down-selected based on available toxicity data and to prioritize chemical diversity for initial screening (Supplementary Table 7). Rett tadpoles or controls were then dosed with lead compounds for 7 days via the culture medium beginning 1–2 days after the onset of symptoms and after collecting baseline swimming behavior data. Animal behaviors (swimming pattern abnormalities and seizure stages) were recorded every other day during drug exposure and washout phases, and the severity of seizure-like events was used as a screening metric.

These studies revealed that the FDA-approved drugs, vorinostat and ivermectin, demonstrated time-dependent reductions in seizure score that reversed after the drugs were washed out (Fig. 2b; Supplementary Fig. 4). However, we focused on vorinostat rather than ivermectin in all subsequent studies due to its higher predicted scoring consistency over time, more favorable reported pharmacokinetics, and its ability, albeit limited, to cross the blood-brain barrier. As a positive control with which to gauge relative efficacy, we also evaluated trofinetide, the synthetic IGF-1 C-terminal tripeptide. Although trofinetide was only moderately effective in Phase III Rett syndrome clinical trials and these results aligned with Phase II trial outcomes¹³, it represents the most relevant clinical-stage benchmark treatment. Importantly, in our Rett tadpole models, both trofinetide and vorinostat reduced seizure-like phenotypes and also increased viability (Fig. 2b). At day 6 of treatment, there was a statistically significant effect by both drugs on reducing swimming score relative to baseline before treatment onset (Fig. 2c). As expected, no significant change was seen in WT and Rett vehicle-treated animals. This head-to-head comparison has supported the possibility that our Xenopus population-level model of Rett can offer predictive value for drug development, which could be further explored using combinatorial treatments in the future.

To understand the protein targets of both drugs, we adapted thermal proteome profiling^68,69,70 that identifies molecules that drugs bind directly to protein targets for use with a whole exothermic animal (Xenopus tadpole) as well as using isolated body segments (head, viscera, tail). While we were able to detect some of the known histone deacetylase (HDAC) targets of vorinostat, we also found that both this drug and ivermectin bind to multiple other proteins related more broadly to acetylation metabolism, beyond histones and HDACs (Supplementary Table 8). These findings suggest a potential mechanism of action for the Rett normalizing responses we observed that involves restoration of normal acetyl-CoA metabolism and post-translational acetylation.

As these observations initially aligned with results of prior studies of HDAC6 acetylation of tubulin in cells of Rett patients^71,72 and microtubule modulation in closely-related CDKL5 deficiency disorder⁷³, we analyzed tadpole tissues for α-tubulin acetylation and observed surprising bidirectional shifts in acetylation patterns depending on the tissue type. For example, a-tubulin was hypoacetylated in neurons in sections of the midbrain (Supplementary Fig. 5a) and hindbrain (Fig. 3f) and increased in GI tract (Fig. 3g) and olfactory multi-ciliated cells (Supplementary Fig. 5b) and following knockdown, resulting in more sparse but dense tangles. Higher resolution imaging revealed that the cilia were longer and misaligned on multiciliated cells compared to controls (Fig. 3a), and both length (Fig. 3b) and orientation (Fig. 3c) were restored by vorinostat treatment. We did not observe a significant effect on cilia morphology in wild-type control tadpoles, which suggests that vorinostat’s impact is disease specific. We then evaluated the impact of ciliary morphology and organization on ciliary function using fluorescent microparticles (Supplementary Fig. 6). The disruption of ciliary beating-induced flow on the surface of skin (Fig. 3d) and olfactory multiciliated cells (Fig. 3e) due to MeCP2 knockdown is reversed by vorinostat treatment, suggesting that ciliary function restoration is linked to ciliary function.

a Olfactory cilia abnormalities due to MeCP2 knockdown were restored by treatment with vorinostat (magenta-yellow colormap of fluorescence intensity relative to z-position, cilia stained for tubulin using mouse anti-α acetylated tubulin and goat anti-mouse Alexa Fluor 594, scale bar = 10 µm), significantly restoring their b length (total number of cilium per condition: N = 43 Rett-Vor and N = 92 Rett-Vehicle, N = 89 Control-Vor, N = 69 Control-Vehicle). and c orientation (total number of multiciliated cells counted per condition: N = 60 Rett-Vor and N = 49 Rett-Vehicle, N = 50 Control-Vehicle). Functional rescue of cilia was also observed in both epithelial, d (N = 3 tadpoles per condition), and olfactory multiciliated cells, e in a microbead clearance assay (total number of tadpoles per condition: N = 4 Rett-Vor and N = 5 Rett-Vehicle, N = 5 Control-Vehicle). Despite the canonical HDACi activity of vorinostat, both hypoacetylated tubulin in the hindbrain (N = 10 Rett-Vor and N = 9 Rett-Vehicle, N = 9 Control-Vehicle tadpoles), f as well as hyperacetylated tubulin in the GI tract (distances from the adventitia epithelium: 26 in 1 µm increments), g was normalized. h Representative immunofluorescent cross-sectional views of the Xenopus gastrointestinal tract stained for ib4 and acetyl-\(\alpha\)-tubulin with Isolectin B₄ FITC Conjugate and mouse anti-alpha acetylated tubulin and goat anti-mouse Alexa Fluor 594, respectively (scale bar = 20 µm; white lines delineate the lumen), and i plot showing that the % of ib4+ cells in the GI tract of tadpoles is increased in Rett tadpoles (N = 4 Rett-Vor and N = 5 Rett-Vehicle, N = 4 Control-Vehicle). Note that the Rett tadpole model has a higher % of ib4+ cells, which is indicative of inflammation and heightened pain response, and that this can be rescued by vorinostat treatment. Error bars indicate standard deviation.

Full size image

Due to the prevalence of GI symptoms in Rett patients and our observation of tubulin acetylation disruption and vorinostat-induced recovery in the GI tract (Fig. 3g), we further explored the role of MeCP2 knockdown in the gut. Staining for isolectin B4-positive (ib4+) cells to assess GI tract inflammation and nociceptive innervation, which is increased in Rett patients^74,75, revealed that tadpoles exhibited significant ib4+ expression due to MeCP2 knockdown (Fig. 3h, i), which vorinostat again reversed.

Vorinostat normalizes the Rett phenotype in a mouse model

Because vorinostat was predicted to reverse the Rett state in our computational studies, we evaluated its efficacy in the MeCP2-deficient male (MeCP2^-/Y) mouse Rett model, previously used in pre-clinical studies for Rett therapeutics, including trofinetide³². We first treated younger animals starting at day 31 post-partum (p31) as was done in previous studies and confirmed that daily intraperitoneal (i.p.) administration of trofinetide (100 mg/kg) results in significant amelioration of multiple disease-related parameters, including diarrhea and motor function, measured by a reduction in the established Bird score¹⁰ that incorporates a broad range of Rett-related CNS and non-CNS symptoms. When vorinostat (50 mg/kg) and trofinetide were similarly administered i.p. on a daily basis over 2 weeks and compared with vehicle control, we found that animals responded favorably to both drugs (Fig. 4a); however, treatment with vorinostat resulted in improved neurological function compared to trofinetide, as determined by measuring elevated plus maze (EPM) performance (Fig. 4b). Both drugs also resulted in enhanced spatial novelty Y-maze performance (Fig. 4c) and an improved diarrhea score (Fig. 4d).

a Bird severity scores measured in MeCP2^-/Y mice treated with vehicle (Rett), trofinetide at 100 mg/kg (Rett-Tro), vorinostat at 50 mg/kg (Rett-Vor (50 mg/kg)) along with wild-type littermate controls treated with vorinostat at 50 mg/kg (Control-Vor) (50 mg/kg) from day 31 to 51 of age (N = 5 mice). Although not indicated in the graph, wild-type littermate controls were treated with vorinostat (N = 3) and trofinetide (N = 3) and did not exhibit signs of severity. b elevated-plus maze (ANOVA and Holm-Šídák test; N = 6 Rett -Vor (50 mg/kg), N = 4 Rett -Vor (100 mg/kg), N = 10 Rett-Vehicle, and N = 9 Rett-Tro, N = 12 Control-Vehicle) and, c Y-novelty maze cognitive tests of MeCP2^-/Y mice (ANOVA and Holm-Šídák test; N = 10 Rett-Vor (50 mg/kg), N = 5 Rett-Vor (100 mg/kg), N = 10 Rett-Vehicle, and N = 7 Rett-Tro) comparing vorinostat treatment efficacy to vehicle and trofinetide, d cumulative diarrhea score during 12 day experimental period using a 0–3 scale (ANOVA and Holm-Šídák test). Representative images of Sholl analysis of microglial arborization in the mouse olfactory bulb, e and the Sholl analysis graph showing representative the branching profiles of microglia analyzed, f in control versus Rett mice treated with or without vorinostat (ANOVA and Holm-Šídák test; N = 5 microglia; Control vs. Rett, P = 0.0004; Control vs. Vor, P < 0.0001). g Levels of indicated cytokines measured in plasma of control mice versus Rett mice with or without treatment with trofinetide (Rett-Tro, 100 mg/kg), vorinostat (Rett-Vor, 50 mg/kg), or vehicle (Rett) for 3 weeks (mean values shown, N = 4 Rett-Vor (50 mg/kg), N = 4 Rett-Vor (100 mg/kg), N = 10 Rett-Vehicle, and N = 9 Rett-Tro, N = 12 Control-Vehicle, scale at right indicates z-score normalized to wild-type vehicle-treated littermates). The observed baseline inflammatory state in Rett mice agreed with published effects of MeCP2 knockdown mice¹⁰⁵. Error bars indicate standard deviation.

Full size image

As brain microglia have been implicated in Rett progression and therapeutic effects^76,77,78 based on their role in neuroinflammation, waste removal, and synaptic pruning, we stained for microglia using and antibody targeting the ionized calcium binding adapter molecule (Iba-1), a microglia-specific marker, in sections of MeCP2^-/y mouse brains. Using Sholl analysis to assess microglia morphology⁷⁹, we found that the MeCP2^-/y microglia exhibited significantly stunted morphology with fewer projection crossings compared to WT littermates, indicative of hyperactivity and a heightened inflammatory state inducing oxidative stress⁸⁰. This was reversed with vorinostat treatment (Fig. 4e, f); however, ramifications in vorinostat treated microglia were higher than in control animals. Together with the reduction of multiple serum inflammatory markers, including TNF-a, IL-4, IL-6, IL-16, CCL2, MIP3b, and CCl1 (Fig. 4g), these findings suggest that MeCP2^-/y mice experience heightened inflammation, observed systemically and in the brain and that treatment with vorinostat can suppress this inflammatory state. Further characterization and classification of the heterogeneity of microglia induced by vorinostat treatment relative to the disease phenotype would be helpful to determine if formulation or further dose optimization is needed.

Oral vorinostat exhibits broad efficacy in Rett mice even after onset of symptoms

Pre-clinical studies for Rett therapies routinely dose MeCP2^-/y mice prior to the onset of symptoms^33,81 due to the rapid onset of symptoms and short lifespan of animals, in contrast to patients who seek clinical care due to developmental deficits which may not be diagnosed for several years before treatment could be initiated⁸². To better mimic this more relevant clinical scenario, male MeCP2^-/y mice were dosed ~1 week after the onset of symptoms. Also, because long term treatment of Rett patients would be best served by oral administration, and due to the fact that our results showed a reduction in diarrhea in mice treated with vorinostat via daily i.p., we tested the efficacy of an oral formulation of vorinostat. Interestingly, when we administered trofinetide i.p. after symptoms developed in this model, it was found to be ineffective based on overall Rett score, similar to the moderate efficacy observed in human clinical trials, whereas oral vorinostat (50 mg/kg) prevented significant worsening of the symptom severity score (Fig. 5a), ameliorated weight gain, and increased performance in EPM (Fig. 5b) and Y mazes (Fig. 5c). Animals treated with oral vorinostat showed complete survival after 3 weeks of treatment, which is comparable to what would be expected for vehicle-treated Rett animals for a study of this duration, whereas only ~60% of trofinetide-treated Rett animals survived (Fig. 5d). To our knowledge, this represents the first effective treatment of MeCP2^-/y mice when treatment is initiated after the onset of symptoms.

a Bird severity scores measured in MeCP2^-/Y mice treated with vehicle (Rett) N = 7, trofinetide at 100 mg/kg i.p. (Rett-Tro) N = 5, oral vorinostat at 50 mg/kg (Rett-Vor) N = 6 from day 36 after onset of symptoms, showing significant suppression of symptoms by vorinostat treatment while trofinetide’s effects were indistinguishable from the vehicle control (ANOVA, ****, P < 0.0001). We also tested wild-type littermate controls treated with oral vorinostat at 50 mg/kg (Control-Vor), N = 5, and similar to previous studies, vorinostat and vehicle treated controls did not exhibit any signs of toxicity or disease severity. b elevated-plus maze and, c Y-novelty maze cognitive tests of MeCP2^-/Y mice comparing vorinostat treatment efficacy to vehicle and trofinetide (ANOVA and Dunnett’s multiple comparison test; N = 7 Control-Veh, N = 6 Rett-Vehicle, N = 5 Rett- Vor). d survival curve of animals in this study. e Cumulative mobility score using a 0–2 scale during the 14-day period and f cumulative diarrhea scored using a 0–3 scale (ANOVA and Holm-Šídák test, N = 7 Control, N = 6 Rett-Vehicle, N = 5 Rett- Vor during a 12-day period). g breathing difficulty evaluated using a 0–1 score (ANOVA and Dunnett’s multiple comparison test; N = 7 Control, N = 6 Rett-Vehicle, N = 5 Rett-Vor). Graphs showing levels of hyperacetylated a-tubulin in multi-ciliated cells in the bronchiolar epithelium of the lung (ANOVA, and Holm-Šídák test, total number of multi-ciliated cells analyzed; N = 14 Rett and Control, N = 16 Rett-Vor) h and in skeletal muscle, i in mice (N = 5 per condition) treated as described in a. j Hematoxylin and eosin (H&E) stained histological sections of colon from the control and drug-treated mice (top row) and immunofluorescent staining for acetylated α-tubulin with rabbit alpha-acetylated tubulin and goat anti-rabbit Alexa Fluor 633 (middle) and βIII-tubulin with chicken anti-βIII- tubulin with goat anti-chicken Alexa Fluor 488 (bottom) in these sections, scale bar = 100 µm. k Graph showing changes in the βIII-tubulin ratio in drug-treated versus control colon tissues (ANOVA and Holm-Šídák test, N = 4 Control, N = 4 Rett-Vehicle, N = 4 Rett- Vor; N = 3 Rett-Tro). Error bars indicate standard deviation.

Full size image

Rett patients and mouse Rett models exhibit significant declines in fine motor and gait control as well as GI tract dysfunction^6,83. Importantly, the deficiencies in mobility, gait, and hindlimb clasping, as well as increased diarrhea we also observed in MeCP2^-/y animals were again rescued by vorinostat treatment (Fig. 5e, f). Breathing also was significantly improved following vorinostat treatment compared to trofinetide and vehicle treatments (Fig. 5g). Given the putative effects of vorinostat on acetylation metabolism and post-translational modification of tubulin and its disruption of ciliary function identified in tadpoles, we also analyzed α-tubulin acetylation in the multi-ciliated cells in mouse lung bronchioles, which revealed an overall increase of hyperacetylated α-tubulin in MeCP2^-/y animals that vorinostat again restored (Fig. 5h). In addition, we explored the potential link of tubulin acetylation with muscle function and found that vorinostat treatment can rescue the disrupted α-tubulin acetylation observed in femoral muscle sections (Fig. 5i). Increased acetylation due to the loss of MeCP2 may play a role in muscle function as hyperacetylation has been shown to increase stiffness and resistance in striatal muscles in vitro⁸⁴.

We continued to evaluate the possible impact of tubulin acetylation by examining the colon, a biomechanically-active organ that is often distended in Rett patients and may play a major role in digestive tract disruption that is a major complaint of these patients clinically⁸³. Initial histological analysis of the colon indicated a greater degree of vacuolization and heightened neutrophil infiltrate in MeCP2-null animals, which vorinostat and trofinetide reversed (Fig. 5j). Staining of these sections for bIII- and acetylated a-tubulin also revealed that acetylated a-tubulin is significantly more colocalized with neuronal staining in the colon in MeCP2-null mice, and that this can be normalized by treatment with either vorinostat or trofinetide (Fig. 5j). However, while vorinostat also normalized the ratio of acetylated a-tubulin to bIII-tubulin, it was unexpectedly worsened by trofinetide (Fig. 5j, k). While the developing GI tract of tadpoles prevented detailed analysis of their enteric neurons, the increased number of ib4+ cells in the Rett tadpole GI tract is consistent with our observation of signs of hyperinflammation in the MeCP2-null mouse colon and with other reports using MeCP2-null mice⁸⁵. The restoration of the number of ib4+ cells and the tubulin acetylation-normalizing effect of vorinostat identified by network analysis in Xenopus clearly translated to this mammalian (mouse) model at the cellular level.

These data show that oral dosing of vorinostat improved classical disease outcome metrics in MeCP2-null male mice and it ameliorated microglial dysfunction and hyperacetylation of tubulin in their GI and respiratory tracts as well as in skeletal muscle. Vorinostat also was recently shown to be effective in treating Fragile X syndrome⁸⁶, which is closely related to Rett, sharing both X-linked mutations and symptoms of autism, although it is caused by a mutation in the FMR1 gene. Interestingly, we identified a loss of gene network connectivity with FMR1 in the Rett-related Xenopus gene network (Fig. 1h), and genes with the maximum variation and fold change differences in humans with Fragile X also occurred in pathway classes involved in oxidative stress and multiple biosynthesis pathways that align with the major regulators we identified by the cross-species MeCP2 mutant network analysis (Supplementary Fig. 7a and Supplementary Tables 9 and 10). Moreover, when we used nemoCAD to predict drugs that might reverse the Fragile X genotype in humans, both the compound classes and gene targets were similar to those predicted for Xenopus with MeCP2 knockdown, including a high number of PPAR receptor agonists, HDAC inhibitors, and cyclooxygenase inhibitors (Supplementary Fig. 7b, c). Thus, these computational results suggest more broadly that developmental therapeutics approaches for neurodevelopmental disorders with similar broad transcriptional effects may benefit from using this form of gene interaction network analysis to identify network-level signatures characteristic of specific diseases and to predict drugs that may reverse those changes.

Current therapeutics for Rett Syndrome focus only on the neurological manifestation of the disease rather than the multiple comorbidities reported in patients⁸⁷. Our AI-enabled computational approach, nemoCAD, tackles this issue by analyzing transcriptomic profiles in an agnostic manner, and by simulating gene to gene and drug to gene networks while looking for mechanistic diversity within a given drug class, which leads to predictions that address the disease as a whole system. These predictions are then screened in our Xenopus Rett model, for whole organism toxicity and efficacy testing, which we created using CRISPR gene engineering technology. These results were then validated by carrying out similar studies in MeCP2 knockdown male mice, MeCP2^-/y. This approach resulted in the discovery that the existing FDA approved drug, vorinostat, is effective at ameliorating key neurological, gastrointestinal, respiratory, inflammatory dysfunction in the Rett Syndrome mouse model. Taken together, our results demonstrate the value of combining AI-enabled, computational network analysis-based drug repurposing algorithms with CRISPR technology to rapidly build in vivo disease models and assess efficacy of drugs with activities of interest predicted computationally.

The Xenopus Rett tadpole model was generated in our laboratory within 3 weeks, and due to the stochastic effects of CRISPR/Cas9 mediated MeCP2 knockdown, our model exhibited a large degree of genetic and phenotypic heterogeneity, which is representative of the disease states reported in the Rett patient population⁸⁸. The Rett tadpole model was analyzed starting at stage 45 dpf when gastrointestinal, respiratory, neural, and skeletal systems are developed⁸⁹, which enabled us to investigate multiple fundamental biological processes simultaneously. At this stage, we observed and characterized CNS and non-CNS related abnormalities, which include bidirectional dysregulation on transcriptional and protein levels, highlighting the multi-functionality of the MeCP2 protein beyond the CNS⁹⁰.

In contrast to existing AI-enabled computational drug prediction frameworks, nemoCAD prioritizes network impact of perturbations rather than a single perturbation signature within a biological system. This allows us to consider response similarity within gene neighborhoods as well as on the level of individual genes for drug prediction without relying on one target. Importantly, other existing analysis frameworks with targeted hits have shown limited success translating predictions for Rett Syndrome beyond in vitro screens, despite incorporating complex knowledge maps that mine published transcriptional datasets^61,91. Furthermore, they were unable to identify our lead compound, vorinostat, as a therapeutic for Rett syndrome, which aligns with our benchmarking results with CMAP⁹². With our agnostic approach, where cell types are not down selected, we generate a ranked list of predicted compounds to improve system-wide perturbations found in Rett Syndrome, which were tested with our Xenopus tadpole disease model to mitigate the risks associated with computational and screening outcomes⁹³.

Our combined AI and experimental approaches enabled us to demonstrate vorinostat’s efficacy from a neuromuscular perspective and prompted the rapid discovery of an underappreciated mechanism: the counterintuitive ability of an HDAC inhibitor to impact acetylation metabolism in a bidirectional manner. This was made possible due to the presence of highly differentiated multi-ciliated cells that are found throughout Xenopus tadpoles, which have been previously used to model kinetics in cerebral spinal fluid, skin epidermis, respiratory system, and gastrointestinal tract^94,95,96,97. In the Rett tadpole model where we observed abnormal swimming behavior and seizures, we also detected widespread dysfunction of α-tubulin acetylation on a morphological and functional level in the respiratory, gastrointestinal, and neurological systems. For example, using immunohistochemical staining, we observed hyperacetylation of tubulin in the enteric neurons of the colon, the multi-ciliated cells in the bronchioles and the skeletal muscle fibers in MeCP2^-/y mice, which corresponded to the respiratory, motor, and gastrointestinal abnormalities noted in a previously described preclinical model^10,85. Importantly, treatment with vorinostat restored these abnormalities in all of these organ systems.

From our predicted list, vorinostat was selected due to its preclinical and clinical profile as a well-tolerated cancer therapeutic⁹⁸. Moreover, vorinostat administration was well tolerated by MeCP2^-/y mice, which experience a severe form of Rett Syndrome as they do not survive past 12 weeks of age¹⁰. When treated before the onset of symptoms in MeCP2^-/y, we observed comparable efficacy with the current FDA treatment, Trofinetide⁸¹. In addition to the neurological improvements previously noted in preclinical studies with Trofinetide^81,99, we observed significant improvement of both gastrointestinal and respiratory dysfunction. When vorinostat was further formulated for oral administration, we continued to demonstrate broad efficacy in MeCP2^-/y mice when treated 1 week after onset of disease in MeCP2^-/y. To our knowledge, this is the first effective treatment to demonstrate broad efficacy dosed after symptoms are presented.

In studies carried out in parallel in the Rett tadpole model, vorinostat restored the hyperacetylation of α-tubulin observed in multiple key organs implicated in Rett. Alpha-tubulin is a ubiquitous cytoskeletal protein that is fundamental to many cellular functions related to microtubule dynamics, including neuronal health, and dysfunction of this protein is associated with several neurodegenerative diseases, such as Parkinson’s¹⁰⁰. Thus, our results provide insight into which organ systems are impacted by this post-translational dysregulation in Rett. Previous studies have explored the impact of hyperacetylation of histones on Rett disease pathology, including histone H3, which capture the subtle bidirectional transcriptional dysregulation that we and other observed with the loss of MeCP2¹⁰¹. Our results suggest that additional studies on the impact that acetylation metabolism have specifically on neuronal, gastrointestinal and respiratory health would prove valuable to further delineate vorinostat’s mechanism of action and the pathogenesis of Rett Syndrome.

While gastrointestinal issues such as diarrhea are reported to be a major side effect with Trofinetide treatment in Rett patients, we did not observe exacerbation of diarrhea when treating MeCP2^-/y mice with Trofinetide. Some have argued that certain ingredients used in the trofinetide formulation may cause the reported side effect in patients¹⁰². Another possible explanation for our findings is that we did not reach a dose high enough to observe this toxicity in mice given that the gastrointestinal issues are reportedly dose-dependent¹⁰³. Nonetheless, Trofinetide treatment did produce a significant increase in α-tubulin hyperacetylation in enteric neurons of the colon. Thus, additional studies may be necessary to explore how this posttranslational modification relates to Trofinetide’s gastrointestinal side effects^102,103 using an escalating dose regimen. But taken together with the modest efficacy of trofinetide in the clinic, our results suggest that conducting pre-clinical studies to restore health after the onset of symptoms as in our Xenopus and mouse models may be a more robust benchmark for efficacy in studies on complex CNS disorders.

In summary, the AI-enabled drug repurposing platform presented here offers Rett and other neurodevelopmental patients the possibility of identifying single small molecule therapeutics within weeks that address multiple fundamental disease-modifying pathways with broad clinical impact. This approach also can enable the discovery and clinical derisking of hidden therapeutic mechanisms for drug development, while providing patients with badly needed treatments. Importantly, vorinostat was recently granted Orphan Drug Designation by the FDA and it is currently moving into human clinical trials for treatment of Rett syndrome based in large part on the data presented in an earlier version of this manuscript.