Transcriptomic Profiling of Zebrafish Mutant for cdkl5 Reveals Dysregulated Gene Expression Associated with Neuronal, Muscle, Visual and Skeletal Development

1. Introduction

Human cyclin-dependent kinase-like 5 (CDLK5) is encoded by a single-copy gene located on the short arm of the X chromosome. It is widely expressed in many tissues throughout the human body. However, it is most abundant in the brain (specifically in forebrain, cerebral cortex, hippocampus, striatum and olfactory bulb), testicles and thymus[1,2,3]. During the embryonic stages, CDKL5 expression is minimal but it quickly increases in postnatal stages and throughout neuronal development and synaptogenesis, indicating that CDKL5 is crucial for early normal brain development and function[3,4,5]. At subcellular level, it is found both in the nucleus and in the cytoplasm, where it plays specific roles involved in important neuronal processes[6].

CDKL5 is an important protein kinase that is associated with several essential functions through the phosphorylation of its targets and interactions or associations with other proteins[7]. Although the extent of CDKL5 exact functions remains to be elucidated, it is known that in the cytoplasm it regulates axon outgrowth through the interaction with shootin1 (SHTN1)[8]. It is also involved in the regulation of the actin cytoskeleton and dendritic arborization through interaction with Ras-related C3 botulinum toxin substrate 1 (Rac1)[9] and in synapse formation via phosphorylation of amphiphysin 1 (AMPH1)[10,11] and the netrin G1 receptor-protein membrane (NGL-1)[5]. CDKL5 also regulates dendritic microtubule dynamics by direct phosphorylation of microtubule-associated protein 1S (MAP1S)[12,13], EB family member 2 (EB2)[13], Rho guanine nucleotide exchange factor 2 (ARHGEF2)[13] and DLG5 (Discs Large MAGUK Scaffold Protein 5)[12]. Furthermore, by direct phosphorylation of CEP131, CDKL5 also plays a role in ciliogenesis[12]. In the nucleus, CDKL5 phosphorylates the DNA methyltransferase 1 (DNMT1)[14] and the methyl CpG-binding protein 2 (MeCP2), thus regulating gene expression through its effect on DNA methylation [15].

Disruption of CDKL5, resulting in decreased activity, leads to CDKL5 deficiency disorder (CDD), a rare and severe condition classified as a developmental epileptic encephalopathy[16,17]. Mutations in the CDKL5 gene are responsible for the development of this disorder, which can manifest with a broad range of symptoms. It is primarily characterized by early-onset, treatment-resistant seizures (which begin within the first three months of life) and severe neurodevelopmental impairments, including deficits in motor, communication, and cognitive development[16,17,18]. Other common phenotypes include hypotonia, cortical visual impairment, microcephaly, subtle dysmorphic features, scoliosis, sleep disturbances and gastrointestinal problems[16,19]. Currently, there is no cure available for CDD affected individuals and since targeted therapies are scarce, treatment only focus on seizure management and in relieving the symptoms[18,19,20].

Despite extensive efforts to elucidate the molecular factors underlying the onset and progression of CDD, its development is still unclear. Several mouse models have been generated to understand how CDKL5 deficiency leads to neurological defects[21,22,23,24,25]. These models recapitulate several phenotypes of the human disorder such as abnormal eye tracking, impaired learning and memory, motor dysfunction and changes in locomotion[21,24]. However, these studies only focus on the neurological phenotypes of the disorder, thus our current knowledge of CDKL5 function in other organs/tissues besides the brain is still limited. The comprehension of affected gene regulation and signaling pathways is crucial for unraveling the mechanisms underlying both normal physiology and disease states. Few studies on transcriptomic analysis in the context of Cdkl5 deficiency have been reported, and those that exist were performed using specific brain regions from Cdkl5 null mice[26]. Previously, we identified zebrafish as a promising tool to investigate CDD. We then characterized and validated the first cdkl5 mutant zebrafish line (sa21938)[27], suggesting its potential as a valuable model for exploring the underlying mechanisms of CDD and its relevance in high-throughput drug screening[28]. This zebrafish mutant line harbors a nonsense mutation that introduces a premature stop codon, leading to either mRNA degradation or the production of a truncated protein, resulting, among other phenotypes, in altered motor behavior, increased susceptibility to seizures, and skeletal defects[28].

The main objectives of this study were to provide new insights into the molecular mechanisms and pathways linked to Cdkl5 functions with a focus on its unexplored roles and to evaluate whether the cdkl5^sa21938 mutant zebrafish line shares molecular similarities with mammalian systems that reproduce the CDD phenotypes. The validation of this model would enable its use to further investigate CDD pathophysiology and contribute to the discovery of new therapeutic molecules. To that end, we investigated the transcriptomic alterations caused by the loss of Cdkl5 function in cdkl5 homozygous zebrafish at five days post-fertilization (dpf) and 35 dpf using high-throughput RNA sequencing approach.

2. Materials and Methods

2.1. Ethics Statement

All animal procedures carry out in the present study were performed in compliance with ARRIVE guidelines (https://arriveguidelines.org) and according to the EU and Portuguese legislation for animal experimentation and welfare (Directives 86/609 CEE and 2010/63/EU; Decreto-Lei 113/2013; Portaria1005/92, 466/95 and 1131/97). This study was approved by the Portuguese Direção-Geral de Alimentação e Veterinária (authorization no. 0421/2021). Qualified operators conducted animal handling and experimentation, prioritizing the minimization of pain, distress, and discomfort.

2.2. Fish Maintenance

Experiments were conducted using the cdkl5^sa21938 mutant zebrafish line (obtained through European Zebrafish Resource Center) and their wild-type (WT) siblings. All adult zebrafish and embryos were maintained in the animal facility in a recirculating water system under controlled temperature (28 °C) and lighting (14h light -10h dark cycle). For the experiments, homozygous cdkl5^sa21938 (cdkl5^-/-) mutant adult fish were incrossed to generate cdkl5^-/- embryos, while WT embryos were obtained by crossing adult WT siblings, which were the offspring of heterozygous cdkl5^sa21938 mutant incrosses. The cdkl5^-/- fish were also crossed with Tg(Hb9:GFP) transgenic reporter fish line[29,30]. The resulting heterozygous zebrafish were incrossed to obtain the homozygous lines, which were confirmed by genotyping.

2.3. RNA Extraction

Total RNA was isolated from embryos and juveniles of homozygous cdkl5^sa21938 and their WT siblings at 5 and 35 days post-fertilization, respectively, using the NZYOL reagent (NZYTech, Lisbon, Portugal) and following the manufacturer’s instructions. Each group comprised five biological replicates, prepared from pools of 50 embryos or 7 juveniles. DNase treatment was then performed using the RNase-Free DNase set (QIAGEN, Hilden, Germany) to eliminate any remaining genomic DNA contamination. To increase purity, RNA samples were subjected to the RNeasy kit (QIAGEN) clean-up protocol. The RNA quantity and purity were assessed using a Nanodrop spectrophotometer (Thermo Fisher Scientific, Waltham, USA.), and its integrity was verified by electrophoresis on a 1% agarose gel.

2.4. RNA Sequencing (RNA-Seq)

The purified RNA was delivered to Novogene company (Cambridge, UK) for library construction and sequencing. Integrity of RNA samples were analyzed on the Agilent 2100 bioanalyzer and only samples with an RNA integrity number (RIN) > 8 were further used for cDNA library preparation. Briefly, mRNA was separated from total RNA using poly-T oligo-attached magnetic beads and subjected to fragmentation using divalent cations under high temperature. Then, the first strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase, followed by the second strand cDNA synthesis using DNA polymerase I and RNase H. The lasting overhangs were transformed into blunt ends through polymerase and exonuclease reactions. To prepare for hybridization, adaptors containing a hairpin loop structure were ligated to the previously 3’ ends adenylated DNA fragments. The AMPure XP system (Beckman Coulter, Beverly, USA) was used to purify the cDNA fragments with 370~420 bp in length that were then digested with USER enzyme. PCR was carried out using Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) primer. The PCR products were purified, and the Agilent Bioanalyzer 2100 system was used to evaluate the library quality. Finally, directional (strand-specific) libraries were sequenced on the Illumina Novaseq platform, generating 150-bp paired-end reads for each sample.

2.5. Bioinformatic Analysis of RNA-Seq Data

For the RNA-Seq data analysis, raw reads in fastq format were first trimmed to obtain clean reads through the elimination of low-quality reads, adapter-containing reads, and ploy-N-containing reads. After filtering, clean reads were aligned to the zebrafish reference genome (Danio rerio GRCz11) using Hisat2 (v2.0.5). For gene expression quantification, featureCounts (v1.5.0-p3) was used to count the number of reads mapped to each gene. Fragments per kilobase per million mapped fragments (FPKM) for each gene were then calculated, considering the gene length and the total number of mapped reads. Genes with FPKM values greater than 1 were considered to be expressed. Analysis of the differential expression between cdkl5 mutant and WT zebrafish was performed using the DESeq2 R package (1.20.0). The resulting p-values were adjusted (padj) according to the method proposed by Benjamini and Hochberg[31] for controlling the false discovery rate. Genes with a p-value < 0.05 and |log2(FoldChange)| > 0 found by DESeq2 were considered as differentially expressed. Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) was performed using the clusterProfiler R package. GO terms with corrected p-value < 0.05 were considered significantly enriched by differential expressed genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to identify significantly enriched pathways associated with DEGs. The clusterProfiler package in R was used to test the statistical enrichment of DEGs in KEGG pathways.

2.6. cDNA Synthesis

Total RNA was reverse-transcribed into cDNA using the Moloney murine leukemia virus (M-MLV) reverse transcriptase (Invitrogen, Carlsbad, USA). Briefly, 1µg of RNA was combined with 0.4 µl of oligo(dT) primer (50 mM), 1 µl of dNTPs (10 mM) and DNase/RNase-free water to a final volume of 12 µl, then incubated at 65°C for 5 min. Next, 4 µl of 5× first strand buffer, 2 µl of DTT and 1 µl of Ribolock RNase inhibitor were added, followed by incubation at 37°C for 2 min. Subsequently, 1 µl of M-MLV was added, and the reverse transcription continued for 50 min at 37°C. Finally, the mixture was incubated at 70°C for 15 min to inactivate the enzyme.

2.7. Quantitative Real-Time PCR (qPCR)

Gene-specific amplifications by qPCR were performed in 20 µl reactions containing 10 µl of SensiFAST SYBR Hi-ROX (Meridian Bioscience, Cincinnati, USA), 0.8 µl of each forward and reverse set of primers (10 µM) and 2 μl of cDNA (diluted 1:10). Gene-specific primers are listed in Supplementary Table S1. The qPCR reactions were carried out on a 96-well plate, using a CFX Connect Real-Time PCR Thermocycler (Biorad, Hercules, USA) and under the following conditions: an initial denaturation step at 95°C for 2 min, followed by 40 cycles of amplification (each cycle is 5 s at 95°C, 20 s at 60°C). After the amplification was completed, a melting curve was generated by gradually heating the sample to 95 °C at a rate of 0.1 °C per second, while continuously capturing the fluorescence. For each sample, at least two technical replicates were performed. Relative gene expression levels were determined applying the ΔΔCt comparative method and were normalized to the β-actin and rps18 reference genes for 5 dpf and 35 dpf samples, respectively (Supplementary Table S1).

2.8. Imaging and Analysis of the Tg(hb9: GFP) Embryos

Embryos with 3 dpf were anesthetized in MS-222, oriented on the lateral side on top of 0.8% agarose and fluorescent images were acquired using the ZEISS AXIO Zoom V16 microscope. Images were captured and processed using ZEN 3.7 software. The deconvolution function (deblurring) was used to improve resolution, using the same parameters in all images. The number of GFP⁺ spinal motor neurons present in three hemisegments after the distal end of the yolk extension were quantified. The length of five motor axons (located immediately after the yolk sac, over the yolk extension) per embryo was measured using the Neuroanatomy/SNT plug-in[32] for ImageJ/Fiji. The average of the five values was calculated for subsequent statistical analysis.

2.9. Statistical Analysis

Statistical analysis was performed using Prism 8 (GraphPad Software). The Kolmogorov-Smirnov test was used to determine the normality of the data. Data with a normal distribution are presented as mean ± SD. One-way analysis of variance (ANOVA) followed by Tukey’s post-test was used to identify significant differences between three groups, while a t-test with Welch’s correction was applied for comparisons between two groups. Data that did not pass the normality test are presented as median with the interquartile region (boxplot). The Kruskal-Wallis test, followed by Dunn’s multiple comparisons test, was used to identify significant differences between three groups, while the Mann-Whitney test was used for comparisons between two groups. Differences were considered statistically significant when p < 0.05.

3. Results

3.1. RNA-seq Data Description and Validation

To identify direct and indirect genes and molecular pathways affected by the potential loss of Cdkl5 function, a high-throughput gene expression analysis through RNA sequencing (RNA-seq) was performed. This analysis compared the transcriptomes of cdkl5^-/- mutant zebrafish and their wild-type (WT) siblings at 5 and 35 days post-fertilization with five biological replicates in each group. The RNA-seq data was deposited in the Gene Expression Omnibus (GEO) repository (accession number: GSE294284).

A summary of the RNA-Seq descriptive data is presented in Table 1.

After sequencing, a total of 466,025,230 and 458,390,836 raw reads were generated from 5 dpf embryos and 35 dpf juveniles’ samples, respectively. The obtained sequencing error rate, the quality score of Q30 and the average GC content for all samples validate the good quality of the sequencing. Subsequent data filtering by removal of low-quality reads, adapter-containing reads, and ploy-N-containing reads from 5 dpf embryos and 35 dpf juvenile, identified the high-quality clean reads, which were aligned against the zebrafish reference genome and mapping rates ranging from 87.40% to 92.41% were obtained (Table 1). The percentage of mapped reads should be higher than 70% when an adequate genome is chosen, and no contamination occurs. Therefore, these outcomes indicate that the obtained sequencing data are accurate and reliable.

3.2. Gene Expression Analysis of cdkl5^-/- Mutant and WT Siblings’ Zebrafish

Expression levels for all identified genes in the cdkl5^-/- and WT samples from 5 and 35 dpf zebrafish were quantified by calculating the FPKM. The distribution of gene expression levels in the analyzed samples is indicated in the boxplot graphs of Figure 1A,B. As shown, the overall range and distribution of the FPKM values were similar and uniform across the five biological replicates of each group at both developmental ages indicating good standardization and reliability of the data.

The correlation of gene expression levels across samples is crucial for ensuring reliability and appropriate sample selection. If this is the case, the square of the Pearson correlation coefficient (R²) should be greater than 0.8. We performed this analysis, and the correlation coefficients (R²) are displayed as heatmaps in Figure 1C and 1D. As shown, the correlation between biological replicates within the cdkl5^-/- and WT groups was high, with R² values greater than 0.911 and 0.930 for 5 and 35 dpf samples, respectively. This indicates that all data within each group were well correlated and that the biological replicates had good reproducibility.

Additionally, principal component analysis (PCA) was performed on the gene expression values (FPKM) of all samples at 5 and 35 dpf, as shown in Figure 1E,F, respectively. The similarity between samples is indicated by the distance between them, i.e. shorter distances correspond to higher similarity. Overall, in both PCA maps for 5 and 35 dpf zebrafish, the gene expression levels of samples within each group were clustered together; that is, the distance between the five biological replicates from each group was relatively small, indicating similar characteristics. In contrast, the samples from the different cdkl5^-/- and WT sibling groups were well segregated, with relatively large separation distances, indicating distinct gene expression profiles between them.

A total of 17,170 and 17,173 genes were found to be expressed in the cdkl5^-/- and WT siblings’ samples of 5 dpf zebrafish, respectively (Figure 1G). Among these genes, 16,641 were co-expressed in the two groups, 532 genes were uniquely expressed in WT siblings and 529 genes were specifically expressed in cdkl5^-/-. At 35 dpf, 17,681 and 17,503 genes were detected to be expressed in the cdkl5^-/- and WT zebrafish, respectively, with 17,045 genes being commonly expressed (Figure 1H), while 458 and 636 genes were exclusively expressed in WT siblings and in cdkl5^-/-, respectively.

Altogether, these results demonstrate that the data obtained from RNA-seq in this study is both reliable and reproducible, thus suitable to be used for further bioinformatics analysis. It also indicates that the transcriptomic profile of zebrafish at different developmental stages is altered when Cdkl5 is defecient.

3.3. Identification of Differentially Expressed Genes Between cdkl5^-/- and WT Siblings’ Zebrafish

To discover candidate genes affected by the lack of a functional Cdkl5, DESeq2 software was used to identify differentially expressed genes (DEGs) between cdkl5^-/- mutant and WT siblings’ zebrafish at both 5 and 35 dpf. A hierarchical clustering analysis based on FPKM data of all DEGs was carried out and heatmaps were generated to observe expression patterns of different genes. An evident separation between the WT siblings and cdkl5^-/- groups were observed, indicating that the cdkl5^-/- mutant fish displayed a unique gene expression pattern, showing many significant differences from the WT fish at both stages of development analyzed (Figure 2A,B).

At 5 dpf, a total of 2,130 genes were found differently expressed between cdkl5^-/- mutant and WT embryos. Among them, 1,178 genes were downregulated, and 952 genes were upregulated in the cdkl5^-/- when compared to their WT siblings (Figure 2C). While at 35 dpf, a total of 4,210 genes were identified as significantly dysregulated, out of which 2,141 genes were downregulated and 2,069 were upregulated in the cdkl5^-/- mutant juveniles in comparison to the WT (Figure 2D). The detailed information of all DEGs at 5 and 35 dpf is described in Supplementary Tables S2 and S3, respectively. Furthermore, to visualize the abundance and distribution of the differentially expressed genes between the two groups, we have generated a volcano plot for each of the developmental stages (Figure 2E,F).

Since the cdkl5 mutant zebrafish line (sa21938) used in this work harbors a nonsense mutation that may lead to activation of the nonsense-mediated mRNA decay (NMD) mechanism, we specifically examined the expression of cdkl5 between WT and cdkl5^-/- groups. Indeed, cdkl5 was one of the genes found to be differentially expressed between the two groups at both stages of development. It was significantly downregulated in both 5 and 35 dpf cdkl5^-/- mutant zebrafish (Supplementary Tables S2 and S3).

The top 40 known DEGs, including the 20 most upregulated and the 20 most downregulated genes in the cdkl5^-/- mutant compared to the WT based on their significance, are listed in Table 2. Our results showed that the most significantly upregulated gene in cdkl5^-/- embryos at 5 dpf was sort1a (sortilin 1a), which encodes a sorting receptor or co-receptor required for the transport of several intracellular proteins from the Golgi apparatus to the cell surface of lysosomes and endosomes[33]. Ligands of sortilin include the neurotensin (NTS), a neuropeptide functioning as neurotransmitter in the brain which is encoded by nts gene[34], that was also one of the most significantly upregulated genes in cdkl5^-/- mutant embryos, suggesting that Cdkl5 may be involved in sortilin/neurotensin pathway. Additionally, guca1c (guanylate cyclase activator 1C) and arr3b (arrestin 3b) genes that have important function in visual phototransduction pathway[35,36] were also upregulated. The most significant downregulated gene in mutant embryos was gamt (guanidinoacetate N-methyltransferase), encoding for an enzyme involved in creatine synthesis and whose deficiency in humans causes a disorder characterized by developmental delay, seizures and hypotonia[37]. Furthermore, other downregulated genes with greater significance include: ctsl.1 and col5a1, involved in collagen catabolic processes and extracellular matrix organization[38,39]; myh7ba, involved in actin filament binding activity and microfilament motor activity being important for muscle contraction[40]; foxj3, involved in the regulation of muscle fiber identity and regeneration of skeletal muscle through the activation of Mef2c transcription[41]; dynlt3 that encodes a member of dynein motor protein complex, involved in the transport of organelles and vesicles toward microtubule[42]; klhdc8a, involved in the induction of primary ciliogenesis[43]; mbd1b, encoding a protein involved in epigenetic regulation and its deficiency is associated with reduced neurogenesis[44]; rs1a, involved in the cellular organization of the retina; lrp2b, encoding a multifunctional receptor involved in the development of many organs such as brain and eye[45]. At 35 dpf, 9 out of the 20 most significantly upregulated genes in the cdkl5^-/- mutant juveniles were related to visual functions. These included crygm2d18, crygm2d15, crygm2d16, crygm2d11, cryba2a, crygm2d13, and crygm2d14 which belong to the crystallin gene family and are involved in lens development and visual perception[46]; rlbp1b which encode a retinoid-binding protein crucial for the correct function of rod and cone photoreceptors, playing several roles in visual cycle[47]; and lim2.4, which encodes an eye lens protein involved in cell junction organization[48]. The two most significantly upregulated genes were olfml3b and kdm5ba. In the brain, olfml3b encodes a microglia-specific protein whose expression increases during neuroinflammation[49]. The gene KDM5B encodes a histone demethylase involved in the regulation of neuronal, bone and muscle development[50]. Regarding the most significantly downregulated genes in 35 dpf cdkl5^-/- juveniles, several were also found to be significantly downregulated in 5 dpf cdkl5^-/- mutant embryos, including gnl3l, nucks1a, pir, slc41a1, ccdc120, ugt1b5, dynlt3, myh7ba and cishb. The top two were opn1lw1 and opn1lw2, which encode two red-sensitive cone opsins critical for phototransduction[51]. Other genes among the top 20 significantly downregulated included myhc4, myom2a, and vgll2b, three genes involved in muscle development[52,53,54]. Additionally, vwa1, which encodes an extracellular matrix protein essential for cartilage structure and development was downregulated; knockout of this gene in zebrafish causes cartilage dysmorphologies[55]. The gene nfasca, which encodes a cell adhesion protein, was also downregulated. This gene plays a key role in neurite outgrowth and fasciculation, as well as in the organization of the axon initial segment and nodes of Ranvier[56]. Mutations in this gene in humans cause a neurodevelopmental disorder with motor dysfunction[57].

3.4. Gene Ontology (GO) Enrichment Analysis of the Differentially Expressed Genes Between cdkl5^-/- and WT Zebrafish

To investigate the potential biological functions associated with the differentially expressed genes found upon cdkl5 mutation in zebrafish, we performed Gene Ontology (GO) enrichment analysis on DEGs between cdkl5^-/- and WT at the two different stages of development analyzed. Gene Ontology is a collection of terms used to describe the functions of genes and classify gene sets. The GO terms are organized into three main categories: biological process (BP), cellular component (CC) and molecular function (MF).

By using the combined list of all 2,130 DEGs found in 5 dpf mutant embryos, GO analysis revealed that the DEGs could be assigned into the three main categories and further classification resulted in the identification of several GO terms. After using a corrected p-value < 0.05, 65 GO terms were considered significantly enriched. Most of them were associated with biological processes (36 terms), followed by cellular component (16 terms) and molecular function (13 terms) (Supplementary Table S4). The ten most enriched GO terms of each category are shown in Figure 3A. Among these categories, the “neuron projection development” GO term belonging to biological process contained the highest number of DEGs (60 genes). In the biological processes, the most significantly enriched GO terms included “myosin filament organization”, “myosin filament assembly” and “striated muscle myosin thick filament assembly”. The most enriched GO terms within the cellular component comprised “actin cytoskeleton”, “myosin complex” and “extracellular matrix”. Regarding the molecular function category, the most significantly enriched GO terms included “extracellular matrix structural constituent”, “actin binding” and “actin filament binding”. GO enrichment analysis of all 4,210 DEGs identified in 35 dpf mutant juveniles showed that 141 GO terms were significantly enriched (padj < 0.05). As in 5 dpf mutant embryos, most of these terms were associated with biological processes (87 terms), followed by cellular components (38 terms) and molecular functions (16 terms) (Supplementary Table S5). The most significantly enriched GO terms in the BP category included “visual perception”, “sensory perception of light stimulus”, “eye development” and “muscle structure development”. Among the CC category, the top significantly enriched GO terms were “contractile fiber”, “sarcomere” and “myofibril”. In the MF category, the most significantly enriched GO terms included “structural constituent of eye lens”, “actin binding” and “protein-containing complex binding” which was also the term with the highest number of DEGs (122 genes) (Figure 3B). Furthermore, several GO terms were similarly significantly enriched in both 5 and 35 dpf cdkl5^-/- zebrafish, mainly related to muscle development, extracellular matrix organization and actin binding. In contrast, the main alterations between 5 and 35 dpf cdkl5^-/- zebrafish were observed in categories related to the eye, which became particularly enriched at the later stage of development.

To have a better comprehension of which biological functions are being activated and repressed, we have also performed GO enrichment analysis of up- and downregulated DEGs separately. Our results showed that at 5 dpf, significantly enriched GO terms (padj < 0.05) were identified for the downregulated DEGs (Supplementary Table S6 and Figure 4A) which were similar to those identified in the combined analysis. For the upregulated DEGs, only the “neuropeptide hormone activity” GO term was found to be significantly enriched using a padj < 0.05 (Supplementary Table S7 and Figure 4B). At 35 dpf, the top significantly enriched GO terms for the downregulated DEGs were again mostly related to muscle development and contraction (skeletal and cardiac), actin cytoskeleton and extracellular matrix (Supplementary Table S8 and Figure 4C). Additionally, GO terms such as “eye development”, “eye morphogenesis”, “central nervous system development” and “skeletal system development” were also found significantly enriched in downregulated DEGs (Supplementary Table S8). For the upregulated DEGs, the significantly enriched GO terms with higher significance were related to visual and sensory perception, lens and eye development and nervous system processes (Supplementary Table S9 and Figure 4D). Other significantly enriched GO terms in upregulated DEGs were associated with transport of diverse ions, activity of ion gated channel, neurotransmitter transport and secretion, and synaptic signaling (Supplementary Table S9).

Overall, these results indicate that Cdkl5 in zebrafish is involved in important functions associated with the regulation of actin cytoskeleton, skeletal and cardiac muscle development, extracellular matrix organization, neuronal development and eye development.

3.5. KEGG Pathway Enrichment Analysis of the Differentially Expressed Genes Between cdkl5^-/- and WT Zebrafish

The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis is used to identify biological pathways that are significantly overrepresented in a set of genes. To further elucidate the pathways affected by the cdkl5 mutation in zebrafish, we conducted a KEGG pathway analysis on DEGs between cdkl5^-/- and WT. Pathways were considered significantly enriched when corrected p-value was less than 0.05.

Separate enrichment analysis of DEGs at 5 dpf revealed three significantly enriched KEGG pathways for the downregulated DEGs, while no significantly enriched pathways were associated with upregulated DEGs. Detailed information on enriched KEGG pathways and involved genes is indicated in Supplementary Table S10. As shown in Figure 5A, downregulated DEGs were associated with “ECM-receptor interaction” (26 DEGs), “focal adhesion” (38 DEGs) and “cardiac muscle contraction” (21 DEGs) pathways, in accordance with what was observed in the GO terms. By enrichment analysis of downregulated DEGs at 35 dpf, eight KEEG pathways were found to be significantly enriched (Figure 5B and Supplementary Table S11). As at 5dpf, the “focal adhesion” (76 DEGs), “ECM-receptor interaction” (37 DEGs) and “cardiac muscle contraction” (28 DEGs) pathways were also enriched in the downregulated DEGs at 35 dpf. Additionally, “regulation of actin cytoskeleton” (55 DEGs), “Wnt signaling pathway” (41 DEGs), “insulin signaling pathway” (36 DEGs), “adherens junction” (34 DEGs) and “ribosome” (31 DEGs) were also identified. For the upregulated DEGs, 16 KEEG pathways were identified as significantly enriched (Figure 5C and Supplementary Table S12). Most of them were related to metabolism, including “drug metabolism-other enzymes”, “retinol metabolism”, “ascorbate and aldarate metabolism”, and “fatty acid metabolism”.

These results further demonstrate that Cdkl5 dysfunction in zebrafish at different stages of development cause downregulation of pathways essential for the interaction between cells and the extracellular matrix, thus affecting tissue organization and development.

3.6. Validation of the RNA-Seq Differential Gene Expression by RT-qPCR

To verify the reliability and accuracy of the differentially expressed genes identified from the RNA-seq results, we randomly selected several DEGs and performed RT-qPCR to assess their expression levels in cdkl5^-/- and WT zebrafish (Figure 6). In 5 dpf embryos, the expression levels of the genes cdkl5, ahsa1a, col1a1b, ndr2, rac3a and nts were analyzed. RT-qPCR results indicated that cdkl5, ahsa1a and col1a1b were downregulated, while ndr2, rac3a and nts were upregulated in cdkl5^-/- embryos, consistent with the RNA-seq results. In 35 dpf juveniles, the expression levels of cdkl5, mmp9, neb, kcna4, mmp13a and olfml3b were analyzed. Our RT-qPCR results showed a significant decrease in the expression of cdkl5, mmp9, neb and mmp13a, while expression of kcna4 and olfml3b increased in cdkl5^-/- compared to their WT siblings, consistent with the RNA-seq results. Altogether, these findings indicate similar expression trends between RNA-seq and RT-qPCR results, therefore confirming the reliability of our RNA-seq data.

3.7. Impact of Cdkl5 Loss on Zebrafish Nervous System

It is widely known that CDKL5 has important roles in the brain and its deficiency in humans causes neurological defects, a main characteristic of CDKL5 deficiency disorder. To further assess the impact of cdkl5 disruption on the nervous system of zebrafish, we specifically investigated the differentially expressed genes associated with neuronal functions between cdkl5^-/- and WT zebrafish. A subset of these DEGs is listed in Table 3.

Several DEGs related to neuronal development were found and many contributed to the enrichment of GO terms at both stages of development. At 5 dpf, GO terms such as “regulation of nervous system development”, “neuron projection development” and “axon development” were significantly enriched in a higher number of downregulated genes compared to upregulated genes. Accordingly, DEGs related to neuronal development were also mainly downregulated at 35 dpf and contributed to the enrichment of GO terms such as “central nervous system development”. Most of these genes play important roles in neuronal morphogenesis, neurite outgrowth, dendritic arborization and axon guidance. For instance, the genes nfasca, ntn1a, nav2b and plxnb2a were downregulated, while stmn1b and pacsin1a were upregulated at both developmental stages. Additionally, at 5 dpf, slit1a, cntf, casp3a, pcdh18b, mycbp2 and zc4h2 were downregulated, whereas draxina, nefmb, nefla, stmn3, stmn4, srcin1b, amigo1, ndr2 and fez1 were upregulated. At 35 dpf, genes that were downregulated included slit2, pax6b, srgap2, prickle1a and prickle2b. Interestingly, neuronal-related CDKL5 known targets or interactors were also identified as differentially expressed. This included the downregulated genes iqgap1 at both developmental stages, arhgef2, dlg5a and rac1b at 5 dpf and smad3a at 35 dpf, as well as the upregulated shtn1 at 35 dpf (Table 3).

Moreover, genes involved in oligodendrocyte development and myelination were found mainly downregulated at both development stages but only contributed to the enrichment of “glial cell development” and “oligodendrocyte development” GO terms at 35 dpf. For instance, the expression of tfeb, notch3, prx, daam2 and plp1b was found to be altered, in both 5 and 35 dpf cdkl5^-/- zebrafish. They were all downregulated, except for plp1b which was upregulated at 5 dpf and downregulated at 35 dpf. In addition, lama2 and actr10 were also downregulated at 35 dpf (Table 3).

Several DEGs associated with synaptic transmission were identified, with these effects being particularly pronounced at 35 dpf, where the upregulated genes significantly contributed to enrichment of many related GO terms including “synaptic signaling”, “neurotransmitter transport”, “neurotransmitter secretion” and “synaptic vesicle exocytosis”. Genes involved in synaptic vesicle exocytosis and neurotransmission release, which is crucial process of neurotransmission, were found mainly upregulated. Among them were syt12, snap25a and stx12l at both developmental stages and syt2a, syt6b, syt11b, rab3ab, rph3ab, rimbp2, unc13a, sv2ca, napbb at 35 dpf (Table 3). Additionally, several DEGs encoding ion channels playing essential role in neuronal transmission[58,59] were dysregulated in cdkl5^-/- at both developmental stages but with a greater number identified at 35 dpf, where upregulated genes also contributed to the enrichment of GO terms related to ion channel activity. These included calcium, sodium, and potassium voltage-gated channels, which are activated by alterations in electrical membrane potential, as well as ligand-gated ion channels like GABA_A, glutamate, glycine and nicotinic acetylcholine receptors, activated by the neurotransmitter binding (Table 3).

Overall, these findings suggest that the cdkl5 disruption impairs neuronal development and synaptic neurotransmission in zebrafish.

3.8. Effects of Cdkl5 Deficiency on Zebrafish Motor Neurons

In our previous work, we showed that cdkl5^-/- mutant zebrafish exhibit an altered swimming performance compared to WT siblings, characterized by shorter distances travelled. Potential causes for the observed phenotype might be related to muscle weakness caused by defects in the neurons innervating the muscle. Therefore, we have investigated the development of motoneurons in the spinal cord by crossing the zebrafish mutants with hb9:GFP transgenic line. The green fluorescence protein (GFP) expression is controlled by hb9 gene promoter, which encodes a transcription factor crucial for differentiation of postmitotic motor neurons, thereby specifically marking primary and secondary motoneurons. Here, we analyzed the number of hb9+ motoneuron cells and the axonal projection length of the caudal primary (CaP) motoneurons at 3 dpf by fluorescence microscopy (Figure 7). Our results showed that the average number of motor neurons (hb9+ cells) per hemisegment is significantly lower in cdkl5^-/- embryos when compared to WT siblings (Figure 7B). Additionally, the axon length of CaP motoneurons in cdkl5^-/- mutants was significantly smaller than that of WT siblings (Figure 7C).

3.9. Impact of cdkl5 Loss on Zebrafish Skeletal System

In our previous phenotypic analyses of the cdkl5^-/- mutant zebrafish, we identified skeletal abnormalities, including craniofacial cartilage defects and impaired bone development[28]. Therefore, we investigated dysregulated genes associated with the skeletal system that might contribute to these phenotypes.

Several skeletal-related genes were found mostly downregulated in cdkl5^-/- and many contributed to the significant enrichment of GO terms such as “skeletal system development”, “cartilage development” and “extracellular matrix” both at 5 and 35 dpf. Furthermore, GO terms such as “bone development”, “ossification” and “regulation of bone mineralization” were only significantly enriched in downregulated genes at 35 dpf. A subset of identified DEGs are listed in Table 4.

Downregulated genes involved in chondrogenesis and cartilage development in cdkl5^-/- zebrafish included its major regulators sox6 at both developmental stages and sox9a at 5 dpf and sox9b at 35dpf, as well as mbtps1, creb3l2 and skia at 5 dpf and chsy1 at 35 dpf. Expression of the well-known osteogenic markers sp7 and runx2a/b was decreased in cdkl5^-/- juveniles, as well as genes involved in mineralization such as tmem119b at 5dpf and tmem119a, ano6 and ptk2bb at 35 dpf. Interestingly, genes with both roles in cartilage and bone development were also identified, including mef2cb at both developmental stages, yap1 and foxl1 at 5 dpf, and runx1, mef2ca, prdm1a/b, panx3, ccn2b, pthlha, nr3c1 and rflna/b at 35 dpf (Table 4).

Dysregulated genes belonging to signaling pathways that are implicated in cartilage development, osteoblast differentiation and bone formation such as Wntless (Wnt), bone morphogenetic protein (BMP) and transforming growth factor beta (TGFβ) signaling pathways were identified (Table 4). In fact, Wnt-related downregulated genes in cdkl5^-/- juveniles contributed to the significant enrichment of the Wnt signaling KEGG pathway at 35 dpf. These genes included wnt7bb, wnt7aa, wnt16, lrp5, lrp6, dkk2, dvl1a, dvl2, fzd1, fzd4, frzb, ctnnb1, ctnnb2, bcl9, tcf7, tle2b/c, nkd1, nkd2a and kremen1 (Table 4).

A wide diversity of genes encoding components of the extracellular matrix (ECM) crucial for the development and maintenance of the cartilage and bone integrity have been identified. Among them, several collagen genes were downregulated in cdkl5^-/- at both developmental stages, including col1a1b and col1a2, which are the main constituents of bone ECM, col5a1, col6a2, col9a1b, col10a1a, col11a1b and col12a1a. The main collagen type in cartilage, col2a1a/b, was also downregulated at 5 dpf. Other key downregulated ECM constituents included the proteoglycans encoded genes dcn at both stages of development, aspn at 5 dpf and hspg2 at 35 dpf; the glycoproteins encoded genes fn1a, fbln2, and emilin1b at 5 dpf and fbln1 and comp at 35 dpf; the thrombospondins encoded genes thbs1, thbs2a, thbs3b and thbs4 at both developmental stages. The expression of spp1 and bglapl, two major bone ECM non-collagenous genes, was reduced at 5 and 35 dpf, respectively, in cdkl5^-/- mutant zebrafish. Furthermore, genes encoding proteins crucial for ECM organization and stabilization through interactions with ECM components were also found to be downregulated, including hapln1 and tgfbi at both developmental stages, matn1 at 5 dpf, and vwa1 and postna/b at 35 dpf (Table 4).

The expression of genes encoding enzymes of the metalloproteinase family responsible for the degradation of ECM proteins, such as matrix metalloproteinases (MPMs) and a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS), was reduced in cdkl5^-/-. These included mmp14a and adamts16 at both developmental stages, mmp15a/b at 5 dpf and mmp2, mmp9, mmp11b, mmp13a, adamts5, adamts7, adamts10 and adamts15b at 35dpf. Moreover, at 5 dpf, genes encoding proteases of the cathepsin family were also downregulated in cdkl5^-/- embryos, including ctsk, ctsl.1 and ctsbb. The expression of the sh3pxd2b gene, which encodes an adaptor protein involved in ECM remodeling, also decreased at 5 dpf (Table 4).

Altogether, these results demonstrate that cdkl5 dysfunction leads to the reduction of genes essential for normal cartilage and bone formation, consistent with the phenotypes observed in cdkl5^-/- mutant zebrafish.

4. Discussion

In humans, pathogenic variants of CDKL5 cause a severe developmental and epileptic encephalopathy, accompanied by multisystemic comorbidities as a consequence of CDKL5 deficiency[17]. However, its molecular mechanisms are still poorly understood, particularly those associated with the less studied phenotypes. To decipher potential cdkl5 molecular mechanisms, RNA-seq was performed in the present study to obtain the whole transcriptomic profile of a cdkl5 mutant zebrafish previously identified as a suitable model for CDD[28,60], both at a more initial and at a later stage of development (5 and 35 dpf). Our findings revealed the dysregulation of several genes in cdkl5^-/- mutant zebrafish along with enriched GO functions and pathways implicated in diverse functions and tissues at both stages of development. The most significant DEGs in both conditions were linked to muscle, neuronal and visual functions, and GO analysis further highlighted the involvement of the DEGs in muscle and neuronal development at both stages and in eye development at 35dpf. These results indicate a Cdkl5 role within these systems in zebrafish, consistent to the associated phenotypes observed in humans. Although the mechanisms underlying the visual dysfunctions in CDD, such as poor eye contact and lack of visual tracking, are not well understood, they are believed to be linked with cerebral visual impairment[61]. Here, we identified the upregulation of genes crucial for visual perception, including crystalline genes encoding structural constituent of eye lens and genes involved in the phototransduction cascade within the retina. Whether this upregulation is due to Cdlkl5 loss or a compensatory mechanism in attempt to counteract potential cerebral visual dysfunctions should be further investigated.

Additionally, GO terms related to ECM organization and actin cytoskeleton regulation were also predominant at both stages. Consistently, KEGG enrichment analysis identified the ECM-receptor interaction and focal adhesion as the most significant downregulated pathways affected by Cdkl5 loss, at both stages of development, and the regulation of actin cytoskeleton pathway was also enriched at 35 dpf. The ECM, focal adhesions and actin cytoskeleton form an interconnected system fundamental for tissue development through the regulation of cell adhesion, migration, proliferation and differentiation, suggesting that their disruption might account for CDD underlying mechanism.

4.1. Deficiency of Cdkl5 in Zebrafish Causes Alterations in Genes Associated with Muscle System

Individuals with CDD exhibit severe impaired gross motor function with abnormal muscle tone (hypotonia)[17]. Previously, we have also demonstrated that the locomotor behavior of this cdkl5 mutant zebrafish was impaired, suggesting muscle weakness caused by disturbances in the neurons controlling muscle system functions and/or in the muscle itself. Accordingly, RNA-seq data showed that loss of Cdkl5 led to the downregulation of many key genes involved in muscle development and differentiation that was consistent throughout both developmental stages. Indeed, functional enrichment analysis highlighted the enrichment of muscle development-related GO terms such as muscle cell differentiation and skeletal muscle tissue development. Therefore, our findings suggest a key role of Cdkl5 for normal skeletal muscle development. Nevertheless, Serrano et al. found no obvious defects in the muscle organization or patterning of this cdkl5 mutant zebrafish at 2 and 6 dpf[60]. Possible reasons for this inconsistency could be the sensitivity of the method used to detect subtle changes or the appearance of defects later in development.

Interestingly, GO and KEGG pathway analysis revealed the enrichment of downregulation genes involved in cardiac muscle development and cardiac muscle contraction, respectively. Although CDKL5 function in the cardiovascular system is poorly investigated, a few studies reported cardiac abnormalities in CDD patients and Cdkl5^+/− female mice, such as prolonged QT interval and arrythmia, [62,63]. Further studies are needed to investigate Cdkl5 function in muscle.

4.2. Cdkl5 Deficiency in Zebrafish Leads to Dysregulation of Genes Involved in Neuronal Functions

Normal brain development requires precise control of neuronal morphogenesis and signaling. Defects in dendrite formation and arborization, axon outgrowth and guidance, and synapse formation result in the dysfunction of neural circuits that are the cause of several neurodevelopmental disorders[64,65,66,67]. A major clinical presentation of CDD is severe neurodevelopmental impairment that affects both cognitive and motor functions. This study revealed that the loss of Cdkl5 in zebrafish led to the dysregulation of several key genes involved in nervous system development at both developmental stages, as well as the contribution of several related-GO enrichment functions. Specifically, we observed both the downregulation and upregulation of genes with key roles in neurite outgrowth and axon/dendritic arborization and in axon guidance. We reported the downregulation of genes such as nfasca, zc4h2, srgap2, prickle1a and prickle2b which are described as positive regulators of neurite outgrowth and branching, as well as mycbp2, which regulates axon extension/guidance and synapse formation[68,69,70,71,72]. Disruption or pathogenic variants in their human orthologs, leading to loss-of-function, have been associated with neurodevelopmental disorders whose phenotypes resemble CDD, including developmental delay, seizures, motor impairment, hypotonia, intellectual disability and facial dysmorphism[70,73,74,75,76,77]. In addition, Cdkl5 loss in zebrafish caused an alteration in the expression of genes encoding proteins known to interact or be phosphorylated by CDKL5 in mammals and playing a role in these processes. These included IQGAP1 and RAC1 which are crucial regulators of actin cytoskeleton that together with CDKL5 form a complex involved in dendritic arborization; DLG5 which plays a role in cell adhesion, cell polarity and transmission of extracellular signals to the cytoskeleton, thereby regulating dendritic spine morphogenesis and synaptogenesis. ARHGEF2, which regulates microtubule and actin cytoskeleton dynamics, as well as focal adhesion, influencing dendritic spine morphology; and SHTN1, which is involved in axon outgrowth[7]. Altogether, our results are consistent with previous studies using other systems such as knockout mice, cultured neurons and human iPSC-derived neurons, where neuronal CDKL5 loss caused defects in axon outgrowth, dendrite arborization and synaptogenesis. Therefore, indicating a conservation in CDKL5 molecular mechanisms between mammals and zebrafish.

Since the reduced swimming behavior of cdkl5^-/- mutant zebrafish might be caused by defects in motoneurons, we have investigated the effect of cdkl5 dysfunction on the motoneurons development using hb9:GFP transgenic line. cdkl5^-/- embryos displayed a reduced number of motoneuron cells and shorter caudal primary axons at 3 dpf, suggesting that loss of Cdkl5 impairs motor neuron development and reduces axonal outgrowth. Accordingly, previous studies showed that CDKL5 knockdown in primary hippocampal neurons and cultured cortical neurons or CDKL5 knockout in hippocampal neurons from mice resulted in the reduction of axons total length[8,9,78]. Moreover, Serrano et al. showed that through the crossing with islet1:EGFP line, this cdkl5 mutant zebrafish presented reduced branching in the middle primary motor neuron at 6 dpf. However, in contrast to our results, no differences in the number of motoneuron cells were observed. This could be due to the distinct markers used to identify motoneurons or the different stages of development that were analyzed. Altogether, defects in motor neurons, along with the defective expression of genes crucial for muscle development, could explain the impaired locomotor behavior in cdkl5^-/- mutant zebrafish.

Genes involved in the oligodendrocyte development and myelination were downregulated in cdkl5^-/- mutant zebrafish from early stages. This glial-cell type is responsible for the formation of myelin sheath wrapping the axons, that is essential for effective neuronal transmission in the central nervous system[79]. Defects in oligodendrocytes development and maintenance impairing myelination have been implicated with several neurological disorders, including leukodystrophies and autism spectrum disorder (ASD) [80,81,82], that is a feature of CDD. These findings suggest that Cdkl5 might play a role in differentiation and maturation of oligodendrocytes and myelination in zebrafish. Interestingly, white matter alterations and abnormal myelination were observed in CDD individuals[16]. Nevertheless, the function of CDKL5 in oligodendrocyte development and myelination process has not yet been investigated.

Precise synaptic transmission is essential for proper brain function. This process begins in the presynaptic neuron, where an action potential opens Ca²⁺ channels. The Ca²⁺ influx triggers neurotransmitters release from the synaptic vesicles (SV) into the synaptic cleft, then interacting with their receptors on the postsynaptic neuron. For this release to happen, the SVs need to be mobilized to the active zone, docked, and primed for membrane fusion[83,84,85]. Abnormalities affecting any of these key modulators of synaptic function have been implicated in the pathogenesis of many neurodevelopmental diseases including epilepsy, intellectual disability and ASD[86]. Previous studies indicated that Cdkl5 loss impaired synaptic transmission in mice[87,88,89,90,91,92]. Consistent with this, our transcriptomic analysis revealed that cdkl5 disruption in zebrafish led to the dysregulation of several genes involved in various steps of the synaptic transmission process, predominantly showing an upregulation pattern. While some of these genes were affected at 5 dpf, most changes occurred later in development, becoming more evident at 35 dpf as the nervous system matures. Genes encoding diverse types of proteins playing crucial roles in the regulation of SV exocytosis and neurotransmission release were overexpressed in cdkl5^-/- mutant zebrafish. Among these were synaptotagmin genes encoding SV membrane proteins, as well as the vesicle docking-regulators rab3ab and rph3ab, the vesicle priming-promoters rimbp2 and unc13a, and the vesicle fusion-contributors snap25a, sv2ca, stx11b and napb. Interestingly, previous study reported a SV2C overexpression in tissue from human epileptic individuals[93]. Moreover, gain-of-function variants in RPH3A and UNC13A have been detected in individuals with phenotypes ranging from epilepsy and developmental delay to ASD[94,95]. These phenotypes are part of the main clinical presentation of individuals with CDKL5 deficiency disorder and an increased susceptibility to seizures was observed in this cdkl5^-/- zebrafish mutant. Therefore, our results suggest that Cdkl5 plays a role in presynaptic function and its deficiency in zebrafish might enhance SV exocytosis, resulting in aberrant neurotransmitter release thus impairing neuronal transmission. So far, few studies have investigated the presynaptic role of CDKL5, and its only known presynaptic target is amphiphysin 1 that is involved in SV endocytosis and recycling[10]. Additionally, it was reported that loss of CDKL5 in primary hippocampal neurons of rats led to slower SV endocytosis, while SV exocytosis remained unaffected[96]. This inconsistency may be due to differences in the specific cell type analyzed, the developmental stage studied, or even species differences.

Ion channels are also known to play a crucial role in synaptic neurotransmission. Their abnormal activity, whether caused by gain- or loss-of-function (GOF; LOF), is associated with numerous diseases classified as channelopathies[97,98]. Here, we identified several DEGs encoding both voltage- and ligand-gated ion channels, most of which were overexpressed in cdkl5^-/- mutant zebrafish. Voltage-gated calcium channels control the Ca²⁺ influx into cells triggered by membrane depolarization. We observed the upregulation of genes such as cacna1aa, cacna1c, cacna1da encoding the calcium channels subunits Cav2.1, Cav1.2 and Cav1.3, whose GOF mutations in humans leads to excessive Ca²⁺ entry and increased neuron excitability causing developmental disorders such as epilepsy or ASD[99,100,101]. Recently, Cav2.3 subunit encoded by CACNA1E was identified as a target of CDKL5 in humans and mice and loss of its phosphorylation leads to channel GOF causing increased calcium influx and increased neuronal excitability[102]. Additionally, we identified the upregulation of voltage-gated sodium channel genes with key roles in the generation and propagation of action potentials, including scn8aa/b, scn1ba/b genes whose mutations in humans are associated with developmental and epileptic encephalopathies (DEE)[103,104,105,106]. Typically, mutations in SCN8A result in gain of function, causing higher sodium currents and channel hyperactivity[104,105,107]. Interestingly, a recent study showed the effectiveness of sodium channel blockers to control seizures in CDD individuals[108]. Furthermore, we have also observed the upregulation of voltage-gated sodium channel genes such as kcna1a, kcnq3, kcnb1, kcnc1b and kcnd1 whose orthologs in humans have been implicated in epilepsy syndromes like DEE by either LOF or GOF mutations[109,110,111,112,113,114]. Accordingly, organoids derived from CDD individuals displayed neuronal hyperexcitability due to voltage-gated channels dysfunction, marked by increased Na⁺ and K⁺ current densities and premature Na⁺ channels opening[115]. This might suggest common molecular mechanisms for the seizure pathophysiology in CDD between fish and mammals. This study also revealed the overexpression of ligand-gated ion channel in cdkl5^-/- zebrafish comprising both glutamate and GABA neurotransmitter receptors, that mediate excitatory (E) and inhibitory (I) synaptic neurotransmission, respectively. The impaired function of these types of receptors leading to disrupted E/I balance is linked with neurodevelopmental disorders[116]. Among identified upregulated genes were the glutamate receptor grin1a and gria3a/b and the GABA receptor gabra1, gabrd and gabrg2 whose both LOF and GOF mutations in human orthologs are known to cause neurodevelopmental phenotypes including epilepsy, cognitive impairment and/or ASD[117,118,119,120,121,122,123]. This study suggests that loss of Cdkl5 in zebrafish enhances glutamatergic and GABAergic synaptic neurotransmission. In accordance, research using CDKL5 mouse models has revealed disturbances in the E/I balance across multiple brain regions. Consistently, overexpression of glutamate receptors has been reported, leading to enhanced glutamatergic neurotransmission[90,91,124]. Selective loss of CDKL5 in forebrain glutamatergic neurons triggered spontaneous seizures in conditional Nex-cKO mice, due to increased excitability in dentate gyrus granule cells (DGCs), accompanied by increased inhibitory activity, which may be a compensatory response to counteract the E/I imbalance[88]. Similar observations were made in CA1 pyramidal neurons, following CDKL5 ablation in forebrain glutamatergic neurons of Cdkl5 Emx1-cKO mice, leading to hippocampal associated learning and memory deficits[92]. In contrast, CDKL5 deficiency in forebrain GABAergic neurons resulted in enhanced glutamatergic transmission and hypercitability in CA1 pyramidal neurons, leading to autistic features[91]. In addition, CDKL5 loss reduced and augmented excitatory and inhibitory transmission, respectively, decreasing the E/I ration in DGCs of Cdkl5^−/y mice thus impairing hippocampal learning and memory[89]. Increased inhibitory GABAergic transmission was also reported in primary visual cortex and perirhinal cortex of Cdkl5^−/y knockout mice, which might underline the CDD visual impairments and novel object recognition memory deficits, respectively[87,125]. Altogether, our findings indicate a CDKL5 role in the regulation of ion channels’ activity, supporting prior studies stating that CDD is in part a channelopathy.

4.3. Loss of Cdkl5 in Zebrafish Disrupts the Expression of Genes Associated with the Skeletal System

Discrete skeletal phenotypes have been observed in CDD individuals, such as dysmorphic facial features and scoliosis[17,19]. Nevertheless, mechanisms underlying their occurrence are still not known, as no investigations have been conducted to specifically address this issue. We previously reported the skeletal abnormalities in the cdkl5^-/- mutant zebrafish[28], and in the present study, we gained insights into the molecular mechanisms underlying these defects.

During development, cartilage is the first skeletal tissue to form, a process that occurs when mesenchymal stem cells (MSCs) differentiate into chondrocytes, which produce an ECM[126,127,128,129]. Bone formation initiates with MSCs and can occur through different forms involving osteoblast's deposition of an ECM and its later mineralization[127,128,129]. In line with impaired cartilage and bone formation, GO enrichment analysis revealed enrichment of downregulated DEGs involved in cartilage development from the earliest developmental stage analyzed, as well as in bone development, ossification and mineralization at the later stage. Loss of Cdkl5 reduced the expression of master regulators of early chondrogenesis like sox9 and sox6, which are crucial for chondrocyte proliferation and differentiation by inducing the expression of cartilaginous ECM components[126]. Consistently, the expression of main cartilage ECM components was impaired including the collagen genes col2a1, col9a1b and col11a1 whose expression is known to be regulated by SOX9[130,131,132], as well as proteoglycans and glycoproteins encoding genes. Abnormal ECM synthesis and composition have been linked to several diseases associated with craniofacial malformation and skeletal defects[133,134]. Overall, the altered craniofacial cartilage structures in the cdkl5^-/- mutant zebrafish might be the result of a decreased expression of crucial regulators of chondrogenesis and cartilage ECM constituents supporting our previous hypothesis that Cdkl5 plays an important role in cartilage formation and development. Moreover, we observed reduced expression of col10a1 and genes encoding ECM- degrading proteases such as mmp13, mmp9 and adamts5 which are well-established markers of terminally differentiated hypertrophic chondrocytes[135,136,137]. This suggests that Cdkl5 loss impairs endochondral ossification, as the expression of these genes in hypertrophic chondrocytes is required for cartilage replacement by bone.

Intramembranous ossification, which occurs through the direct differentiation of MSCs into mature osteoblasts capable of producing bone ECM[129], also appears to be regulated by Cdkl5 in zebrafish. Loss of Cdkl5 resulted in decreased expression of key osteogenic markers, including runx2 and sp7, which encode transcription factors crucial for early differentiation from MSCs to pre-osteoblast, as well as major bone ECM-encoding genes such as the early osteoblast marker col1a1b, the middle to late marker spp1 and the late marker bglapl[138]. Signaling pathways involved in skeletal development were affected in cdkl5^-/- mutant zebrafish, particularly the Wnt signaling pathway, which was significantly downregulated at 35 dpf, as determined by KEGG analysis. Studies in human and animal models have demonstrated that Wnt signaling induces bone formation, and its dysregulation is associated with bone disorders such as osteoporosis[139,140,141,142]. These findings indicate that the impaired or delayed mineralization phenotype of cdkl5^-/- mutant zebrafish results from defective osteogenic differentiation, suggesting that Cdkl5 plays an important role in bone development.

5. Conclusion

In summary, this study contributes to gain insights into the Cdkl5 molecular mechanisms through the identification of genes, functions and pathways affected by its disruption. Overall, our findings highlight the key role of Cdkl5 in neuronal development and synaptic neurotransmission, suggesting a common Cdkl5 function between zebrafish and mammalian models. Moreover, we demonstrate that cdkl5 disruption causes motor neuron defects that in combination with abnormal muscle development and structure may contribute to the observed gross motor impairments. Finally, this study provides novel evidence of Cdkl5 role in cartilage and bone development, potentially underlying the skeletal defects observed in mutant zebrafish and individuals with CDD. Altogether, these results emphasize the cdkl5^-/- mutant zebrafish as a valuable model for CDD and highlight its potential for screening novel drugs to treat the disorder’s symptoms.