The Zebrafish Core Facility at the Luxembourg Centre for Systems Biomedicine (LCSB) is registered as an authorized breeder, supplier, and user of zebrafish with Grand-Ducal Decree of 20 January 2016 and 26 January 2023. All practices involving zebrafish were performed following European laws, guidelines, and policies for animal experimentation, housing, and care (European Directive 2010/63/EU on the protection of animals used for scientific purposes) and following the principles of the Three Rs. Moreover, we carefully comply with the ARRIVE guidelines. Authorization number LUPA 2020/21 allowed the generation and characterization of the pex1 mutant line, and authorization number LUPA 2017/04 allowed the performance of fin biopsies for genotyping purposes.

Wild-type (AB) and pex1 mutant zebrafish lines were kept in the LCSB Aquatic Platform Core Facility (RRID:SCR_025429) according to standard protocols (Westerfield, 2000; Aleström et al., 2020). Zebrafish embryos were obtained by natural spawning and reared in 0.3X Danieau’s solution [17 mM NaCl, 2 mM KCl, 0.12 mM MgSO₄, 1.8 mM Ca(NO₃)₂, 1.5 mM HEPES pH 7.5, and 1.2 μM methylene blue] in 14 h light–10 h dark condition at 28 °C (± 0.5). Adult fish were maintained in 3.5 liter tanks at a stock density of 4–6 fish/liter with the following parameters: water temperature, 28.5 °C (± 0.5); light:dark cycle, 14 h:10 h; pH 7.5; conductivity, 700–900 μS/cm. Fish were fed three to five times per day, depending on age, with granular food (SDS diets) and live feed (Artemia salina and/or L-type rotifers). For live imaging, larvae were first anesthetized with buffered 0.008% tricaine methane sulfonate (MS222, Sigma-Aldrich) and then embedded in 2.5% methylcellulose. Euthanasia was performed by an approved method (i.e., overdose of anesthetic) or whenever required in larvae (for metabolomics) by rapid cooling at temperatures of 2–4 °C (approval numbers LUPA 2017/03 and LUPA 2022/06).

Multiple sequence alignment was performed with the MUSCLE (RRID:SCR_011812) software from EMBL-EBI using the following peroxisome biogenesis factor 1 protein sequences: Danio rerio (NP_001164306.1), Homo sapiens (NP_000457.1), Mus musculus (NP_001280735.1), and Drosophila melanogaster (NP_652016.1). Jalview software (RRID:SCR_006459, v2.11.4.1) was used to prepare the alignments for visualization.

Candidate guide RNAs (gRNAs) targeting the zebrafish pex1 gene were identified and selected using the online tool CHOPCHOP (RRID:SCR_015723)¹. The generation and purification of gRNAs were performed using the Ambion T7 MEGASCRIPT kit and Qiagen PCR Cleanup kit (#28106) as described in Gagnon et al. (2014). Gene-specific oligonucleotides containing the T7 promoter (5′-TAATACGACTCACTATA-3′), the 20 bp target site (5′-AGATGGCTGACAGACCTAGG-3′), and a 23-bp complementary region (5′-GTTTTAGAGCTAGAAATAGCAAG-3′) were annealed to a constant oligo (AAAAGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACG GACTAGCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC). Microinjections using wild-type (AB) zebrafish and subsequent generation of the mutant line were performed following the same pipeline as described in Heins-Marroquin et al. (2024). 200 ng of pex1 gRNA (0.5 μl), 200 ng of slc45a sgRNA (0.5 μl), 160 ng of Cas9 Nuclease NLS from Streptococcus pyogenes (1 μl, NEB), and 0.2 μl phenol red were mixed in an Eppendorf tube and incubated for 5 min at room temperature. The mixture was placed on ice, and 1 nl was injected into one-cell-stage embryos. For confirmation of the mutation, the CRISPR target site was PCR-amplified from genomic DNA of F1 larvae using NEB Phusion polymerase (NEB, #M0530S) and the following primer pair: zf-5-2-s (5′-CAAAGGAGAAAATTGCTGTTCC-3′) and zf-5-2-as (5′-CTCCGTGTCCTTCTCTTTATGG-3′). The mutant band was cloned into the Promega pGEM-T Easy vector (Promega #A1360) and sequenced in both directions using M13 primers. Sequences were aligned to the pex1 reference gene using the SnapGene (v6.2.2) software. F1 fish carrying the 17 bp deletion were incrossed to obtain stable F2 generations. Genomic DNA from F2, F3, and F4 generation larvae and adult fish was analyzed in the same way as described for the F1 generation.

The DNA from 4 and 5 dpf larvae was extracted following the protocol described in Zhang et al. (2020). Single larvae were placed in wells containing 40 μl extraction buffer (EB; 30 mM Tris-HCl pH 8, 25 μg/ml Proteinase K) in a 96-well plate, and shaken at 220 rpm for 40 min at 37 °C. The EB (≈40 μl) was transferred into a PCR 96-well plate, mixed with 15 μl lysis buffer (10 mM Tris-HCl pH 8.0, 50 mM KCl, 0.3% Tween 20, and 1 mM EDTA), and incubated at 98 °C for 5 min. 7 μl of this mixture was used for PCR amplification using the DreamTaqTM Green PCR master mix (K1082, Thermo Fisher) and the primer pair (zf-5-2-s/zf-5-2-as). Amplicons were analyzed by electrophoresis in 3.5% agarose gels and larvae were grouped according to their genotype.

To monitor survival, a pool of 50–57 larvae at 7 dpf were placed in 1 liter plastic cages (Tecniplast) containing 300 to 500 ml system water and fed daily ad libitum with rotifers. After 6 days, artemia was included in the diet. Daily visual assessment was performed, and dead or severely cachectic larvae were recorded and removed from the cage (and euthanized if needed).

Adult zebrafish were euthanized in 0.04% MS222 (buffered to pH∼7.5) at 4, 7, or 10 months of age. The total length of each fish was measured from the snout to the border of the base of the tail (without caudal fin) while the fish was lying on its side. Fish were then photographed, carefully placed on a dissection mat, and organs were collected using dissecting forceps and surgical scissors. Organs were briefly rinsed in chilled 1X DPBS (no calcium, no magnesium, Gibco, #14190250), and either fixed or snap-frozen in 2 ml Precellys tubes in liquid nitrogen. The tubes containing the snap-frozen organs were weighed and stored at −80 °C until further processing. Fixation was done in either 2% or 4% paraformaldehyde (PFA, Sigma-Aldrich, #252549) overnight at 4 °C.

Locomotor activity of 7, 11, 13, and 15 dpf larvae was recorded using a DanioVision system equipped with a tapping device and a temperature control unit (Noldus Information Technology). The EthoVision XT 11.5 software (RRID:SCR_000441, Noldus Information Technology) was used to analyze the movement of individual larvae. Motion values were expressed as mean swimming velocity (mm/s). Recording protocols used light, dark-light-dark programs, and acoustic-vibrational tapping as stimuli. One day prior to the recording, zebrafish larvae were carefully transferred into each well of a 24-well plate, with 1 ml Danieau’s medium per well. For locomotor recording, each plate was placed in the DanioVision system acclimated at 28 °C and the larvae were allowed to habituate to the new environment for 15 min. After this period, locomotor activity was continuously tracked for 1 h under standard light conditions. Next, the tapping program started with ten tap stimuli (intensity setting: 6) at intervals of 30 s. After a recovery period of 5 min, the light-dark program was initiated with alternating intervals of 15 min. All recordings were performed between 2 and 6 pm. Tracking experiments were performed in independent triplicates (with at least 24 larvae per genotype per replicate). Technical outliers (i.e., larvae not detected by the system) were removed from the quantification. Statistical analysis was performed using the GraphPad Prism analysis software (v10.2.1).

11 and 13 dpf larvae were euthanized in 0.04% MS222 (buffered to pH∼7.5) and fixed in 4% paraformaldehyde overnight at 4 °C. The next day, fixed larvae were washed in 1X DPBS (no calcium, no magnesium, Gibco) and stained in 0.3% Oil-Red O (Sigma-Aldrich, #O1391) as described in Takashima et al. (2021).

The eyes and the liver of adult zebrafish were harvested and fixed as described above. Prior to usage, the adult eyes and livers were washed in 1X PBS (Gibco) and pre-embedded in 15% sucrose solution (in 1X PBS) for 1 h at room temperature for cryoprotection, and then in 30% sucrose solution (in 1X PBS) overnight at 4 °C. The organs were then embedded and frozen in Optimal Cutting Temperature (OCT) (VWR, Radnor, PA, USA) compound in specific cryomolds containing the orientation information (VWR, #720-0821).

Retinal cryo-sections (5 μm) (Leica CM1800 Cryostat) were blocked (0.1% Triton, 10% BSA in 1X PBS) for 1 h at room temperature, washed in 1X PBS, incubated at 4 °C overnight with primary antibody in incubation buffer (0.1% Triton in 1X PBS), washed, incubated at room temperature (RT) for 90 min with secondary antibody, washed, incubated 15 min with DAPI and washed again. Primary and secondary antibodies are described in Supplementary Data 1. Coverslips were mounted using ProLong Gold antifade reagent (Invitrogen, Carlsbad, CA, USA) and retinas were visualized using a Zeiss LSM780 confocal microscope and Zeiss Zen Black software version 2012. Retinal layer measurements were performed using ImageJ (RRID:SCR_003070). For each microscopy image, the thickness of the retinal layer was quantified at 20× magnification. Images were analyzed using ImageJ software. The thickness of the outer nuclear layer (ONL) and inner nuclear layer (INL) was measured at three equidistant points along the retinal section, spanning from the dorsal to the ventral edge. At each point, a straight line was drawn perpendicular to the retinal layers, and the thickness was recorded using the “Measure” function. The mean and standard deviation of ONL and INL thickness were calculated for each group, and data were presented as the average thickness per layer across all measured points. Glutamine synthetase (GS) intensity and retinal layer thickness based on GS-labeled cells were quantified using a MATLAB script (RRID:SCR_001622) designed for image analysis (version R2024a). The script measures the fluorescence intensity and calculates the length of glutamine synthetase-positive regions within the images. Similarly, cone density was quantified as the number of cones per unit area using the same MATLAB script. These automated analyses ensured consistent and reproducible quantifications across all samples.

Fixed livers (4% PFA at 4 °C overnight) were dehydrated with increasing concentration of ethanol (30%–100%) for tissue preservation. Livers were placed in cassettes (Epredia, Fisher Scientific, #10-009-57) and fitted in an automated paraffinization machine (Leica ASP300S) for dehydration, clearing, and inclusion of paraffin into the organ. After embedding, the paraffin block was sectioned with a microtome (Leica RM2155) into 10 μm slices. The latter were dipped in water at 42 °C and placed on cover slides to adhere. Slices were deparaffinized with xylene and decreasing ethanol concentrations (100%–30%) for 5–10 min each. The slides were quickly washed and then incubated in Mayer’s Haematoxylin (Rowley Biochemical, #L-756-1A) for 15 min, followed by rinsing in warm water and washing for a minute. After that, slides were stained in eosin for one minute and dehydrated with increasing concentrations of ethanol and, lastly, xylene. The slices were mounted with a coverslip using Permount (Fisher Scientific, #SP15-100). Livers were visualized using a Leica DMI600 microscope with DFC345FX camera and LASX software (Richmond Hill, Canada). Quantification of peroxisomes in zebrafish liver sections was performed using ImageJ. Images were converted to 8-bit grayscale and subjected to Gaussian blur to reduce noise. Background was subtracted using the rolling ball algorithm (∼20 pixels). Manual thresholding was applied to highlight peroxisomal puncta. The image was then converted to a binary mask, followed by erosion and watershed separation to distinguish closely apposed peroxisomes. Peroxisomes were quantified using the Analyze Particles tool, with size set to 10–200 pixels and circularity 0.10–1.00. The resulting count of peroxisomes was normalized to the number of DAPI-stained nuclei within each image. All quantitative analyses were performed on Average Intensity Projections to reflect peroxisomal content across the entire tissue sample. Data from multiple fields and biological replicates were compiled and analyzed in GraphPad Prism.

Larvae at 11 dpf were euthanized by rapid cooling (2–4 °C) and pools of 10 larvae collected in 0.5 ml Precellys tubes. Total RNA was extracted using Trizol as previously described in Heins-Marroquin et al. (2024) and RNA sequencing was performed by the Beijing Genomics Institute (BGI) using DNBSEQ™ sequencing technology platforms (Drmanac et al., 2010). The raw RNA-seq fastq files were quality-checked using the FastQC (RRID:SCR_014583) software (Babraham Institute Bioinformatics Group. FastQC, v0.12.0, 2023), and the data was pre-processed using the software package Rsubread (R Core Team, 2018; Liao et al., 2019). Gene-level differential expression analysis was conducted in the R statistical programming software (R Core Team, 2018) using the “DESeq2” (RRID:SCR_015687) package (Love et al., 2014) and filtering out genes with low expression counts using the filterByExpr-function with default parameters from the package edgeR (Robinson et al., 2010). P-value significance scores for differential expression were computed using the Wald test, and were adjusted for multiple hypothesis testing according to the Benjamini and Hochberg method (Benjamini and Hochberg, 1995) Pathway enrichment analyses were conducted with the clusterProfiler R package (RRID:SCR_016884, v3.16.0) (Yu et al., 2012) using Fisher’s Exact Test and the gene-level differential expression analysis results obtained with DESeq2 as input. The process network enrichment was visualized in a dot plot using React JS (v19). The pathway over-representation analysis statistics, including false-discovery rate (FDR) scores according to the method by Benjamini and Hochberg (Benjamini and Hochberg, 1995), were determined for gene sets from the Gene Ontology database using gene annotations specific to zebrafish, obtained from the Zebrafish Information Network (ZFIN) database (Sprague et al., 2003).

Ten-fold diluted cDNA samples (2 μl) were added to a mixture (10 μl total volume) containing iQ SYBR Green Supermix (5 μl; Bio-Rad) and the primers of interest (Supplementary Data 1) at a final concentration of 0.25 μM. Reactions were run in a 384-well plate using a LightCycler^® 480 system (Roche). Expression levels were calculated relative to WT samples using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Rpl13α was used as a reference gene for normalization.

The zebrafish genome encodes only one ortholog candidate of the human PEX1 protein (NP_000457.1). The putative zebrafish Pex1 protein (NP_001164306.1) comprises 1237 amino acids and presents 48% and 64% sequence identity and similarity to its human counterpart, respectively. Multiple sequence alignment of peroxisome biogenesis factor 1 proteins from human, mouse, zebrafish, and Drosophila confirmed high conservation of the two N-terminal domains (NTDs) and the two catalytic AAA-ATPase domains, supporting evolutionary conservation of protein function (Supplementary Figure 1). We generated a zebrafish pex1-deficient line using CRISPR/Cas9, targeting the fifth and longest exon of the gene (Figure 1A). This line (institutional code pex1^lux4, following Zebrafish Information Network guidelines) carries a 17 bp deletion in the pex1 gene, causing a frameshift and a premature stop codon in both genomic DNA and mRNA (Figures 1B, C). The predicted truncated protein (305 amino acids) lacks the C-terminal AAA-ATPase domains essential for Pex1 function (Figure 1D). Because antibodies for zebrafish Pex1 are unavailable, and we found that antibodies against human PEX1 did not cross-react, the absence of Pex1 protein could not be verified by immunoblotting. However, proteomic analysis confirmed the complete absence of the Pex1 protein in the pex1^–/– samples, while other Pex proteins were detected, validating the successful knockout of Pex1 in our zebrafish model (Figure 1E and Supplementary Data 1).

To characterize this newly generated line and to determine whether pex1 deficiency is lethal in zebrafish, F2 heterozygous (pex1^+/–) fish were incrossed to obtain F3 offspring, yielding an expected 25% homozygous mutant (pex1^–/–) progeny. A comparable survival rate was observed between pex1^–/–, pex1^+/– and WT siblings in the F3 offspring until 15 dpf, suggesting that pex1 is not essential at early developmental stages (Figure 2A). We suspect, however, that Pex1 function plays a critical role during juvenile stages, as pex1^–/– fish represented only 7%–10% of the population (instead of 25%) by 90 dpf (Figure 2B). Notably, these pex1 knockout adult “survivors” did not show any overt phenotype. However, their lifespan was significantly reduced, and spontaneous death was observed between 6 and 10 months. A humane endpoint was applied to all fish (WT and pex1^–/–) once spontaneous deaths occurred in the pex1^–/– population for three consecutive days. In addition, we suspect that pex1^–/– animals are infertile, as repeated crossing attempts failed to produce any viable offspring. The pex1^–/– females seemed to have functional ovaries producing eggs, which were, however, not released via natural spawning, eventually leading to gonad degeneration (Figure 2C). No significant differences were observed in body length between 7-month-old pex1^–/– and WT fish (Supplementary Figure 2A). Interestingly, both pex1^–/– males and females had developed fatty livers by 4 months of age. These livers typically presented a dense and yellowish appearance and, in some cases, were markedly enlarged, distinguishing them clearly from the reddish livers of WT siblings (Figure 2C). Histological analysis of liver sections confirmed massive macrosteatosis in 7-month-old pex1^–/– fish (Figure 2D).

To assess peroxisome biogenesis in our pex1 mutant model, we immunostained for two complementary peroxisome markers in 7-month-old zebrafish livers: catalase, a matrix enzyme that relies on the PEX5-mediated PTS1 import pathway, and PMP70 (ABCD3), an integral membrane protein that is inserted into peroxisomes via the PEX19/PEX3/PEX16 route and therefore does not require PEX5 (Fang et al., 2004; Okumoto et al., 2020). As expected, catalase staining in WT hepatocytes appeared as fine, punctate foci typical of peroxisomes. In pex1^–/– hepatocytes, however, the catalase signal was more diffuse throughout the cytosol, indicating a failure of matrix-protein import (Figures 3A, B). PMP70 still formed puncta in mutant tissue, confirming that membrane proteins can be incorporated independently of Pex5, yet the number of these puncta was markedly reduced compared with WT. This reduction points to a lower peroxisome abundance rather than a complete loss of the organelle.

To probe peroxisomal metabolism in the pex1 mutant model, we measured fatty acids in liver and brain extracts of ≈7-month-old fish by LC-MS. Heterozygous and WT siblings showed similar levels across all analytes, confirming that a single wild-type pex1 allele sustains normal peroxisomal function. In contrast, pex1^–/– tissues accumulated high levels of phytanic acid, whose α-oxidation is exclusively peroxisomal, and polyunsaturated VLCFAs (C24 and C26) and ultra-long-chain fatty acids (ULCFAs, ≥ C27), reaching log₂ fold changes > 8. These lipid signatures demonstrate impaired peroxisomal α- and β-oxidation, corroborating the loss of functional peroxisomes in pex1^–/– fish (Figure 3C; Supplementary Figures 2B–D, 3 and Supplementary Data 1).

Interestingly, unlike in ZSD patients (Takashima et al., 2017), the brains of pex1^–/– zebrafish showed significantly lower levels of saturated C22 and C24 (Supplementary Figure 2B), and overall saturated FA levels were consistently decreased in pex1^–/– brain and liver samples, in stark contrast to the increased PUFA levels (Supplementary Figures 2B–D, 3). A notable exception among the PUFAs was docosahexaenoic acid (DHA, C22:6n-3), whose levels were decreased in the brain and in the liver of pex1^–/– fish compared to WT (Supplementary Figures 2B, 3). Although this difference did not reach statistical significance, the significantly increased levels of C24:6 (peroxisomal DHA precursor) are consistent with metabolically impaired peroxisomes. The more moderate than expected decrease in DHA levels [Zellweger syndrome patients show strongly decreased tissue DHA levels (Martínez, 1990)] may be explained by an alternative pathway for DHA biosynthesis reported in teleost fish, known as the “Δ4 pathway” (Oboh et al., 2017), or by dietary factors, potentially compensating, at least in part, for deficient peroxisomal DHA biosynthesis. To investigate the latter possibility, we analyzed the fatty acid content of the fish food (Supplementary Figure 4). The SDS100 dry food formulated for feeding in early developmental stages (5 dpf to max. 30 dpf) showed in general higher levels of FAs and DHA compared to SDS400 (dry food for adult maintenance) and live feed. Among the live feed, rotifers, intended for feeding young larvae, contained ten times more DHA than Artemia nauplii, the live feed for juvenile larvae and adult fish. Overall, this fish feed is relatively rich in DHA. For instance, SDS400 and rotifers contain about 10 mg of DHA per gram (≈1% w/w), a concentration in the same order of magnitude as that found in the most DHA-rich human foods, i.e., cold-water oily fish such as mackerel, herring, and salmon (Mozaffarian and Rimm, 2006). This suggests that the zebrafish nutrition may partially compensate for the lack of DHA biosynthesis in the pex1^–/– larvae and adults. Except for the more moderate than expected decrease in DHA levels, our observations in adult pex1^–/– zebrafish suggest that Pex1 plays an important role for proper peroxisomal function in zebrafish.

To identify the earliest phenotypes in our pex1^–/– zebrafish model, homozygous mutant larvae were visually evaluated under brightfield light microscopy. They showed completely normal development until 7 dpf, with no significant difference in body length when compared to WT siblings. However, 96% of pex1^–/– larvae at 11 dpf showed darkening of the liver, and hepatic steatosis was confirmed using oil red (ORO) staining at 13 dpf (Figures 4A, B). Accumulation of polyunsaturated VLCFAs (C24–C30) and methyl-branched fatty acids was confirmed in 11 dpf pex1^–/– larvae using LC-MS (Figures 4C, D and Supplementary Data 1). The most prominently elevated FAs were C26:6, C26:7, C28:8, phytanic acid, and pristanic acid, which showed a log₂ fold change higher than 2 in pex1^–/– larvae compared to the WT siblings.

To further characterize alterations in lipid homeostasis, we performed a targeted lipidomics analysis of lipid species with C14–C22 acyl chains (Supplementary Data 1). Our analysis revealed that, among the lipids analyzed, triacylglycerides (TGs) were the main altered lipid class in pex1^–/– larvae (136.4 ± 53.40 nmol/mg DNA), with about 2 times higher total TG levels than in WT larvae (70.2 ± 20.33 nmol/mg DNA) (Figure 4E). There was a significant increase in polyunsaturated C20 and C22 FAs, including C22:6, and a lower abundance of shorter, saturated and mono-unsaturated FAs such as C16:0 and C18:1, in the composition of TGs, in pex1^–/– larvae compared to WT siblings (Supplementary Figures 5A, B).

Besides TGs, there was a slight but significant increase in hexosylceramides, dihydroceramides, phosphatidylcholine (PC), and phosphatidylcholine plasmalogen (PC-PL) in the pex1^–/– larvae; however, no major alterations in FA composition were observed in these lipid classes (Figure 4E). Interestingly, a consistent decrease in phosphatidylethanolamine (PE) and phosphatidylethanolamine plasmalogen (PE-PL) was observed in the pex1^–/– larvae, with only 1-alkyl,2-acylphosphatidylethanolamine being significantly lower (Figure 4E). For the remaining analyzed lipids, including diacylglycerides (DG), phosphatidylinositols (PI), phosphatidylserines (PS), phosphatidylglycerols (PG), cholesterol esters (CE), ceramides (CER), lactosyl-ceramides (Lac-Cer), lysophosphatidylethanolamines (LPE), and lysophosphatidylcholines (LPC), no significant changes were found in the pex1^–/– larvae versus WT siblings (Supplementary Figure 5C).

We next performed bulk RNA sequencing on pex1^–/– larvae and WT siblings at 11 dpf. Partial least squares discriminant analysis (PLS-DA) separated the samples by genotype along the first component, indicating that the pex1 mutation produces a distinct transcriptomic profile (Figure 5A). Differential gene expression analysis revealed 566 differentially expressed genes (DEGs) (p-adjust < 0.01), with 415 upregulated and 151 downregulated DEGs in pex1^–/– versus WT larvae (Figures 5B, C and Supplementary Data 1). Six DEGs (rida, arf4b, hpxa, isg15, scd, and hspa5) were validated by qPCR, corroborating the transcriptomic findings (Figure 5D). Among the DEGs, 12 downregulated and 32 upregulated genes displayed a log₂ fold change of more than 1.5. Among these, arf4b, isg15, scd, ppp1r15b, fosl1a, and plin2 were the most upregulated, while the liver-specific genes hpxa, ahsg1, and cyp7a1 were the most downregulated genes in the pex1^–/– larvae. These data implicate the liver as an early, and possibly initiating, site of pathology in this model.

To identify altered pathways potentially relevant for Zellweger syndrome, DEGs were mapped to their human orthologs according to the HUGO Gene Nomenclature Committee (HGNC). Pathway enrichment analysis was then performed using the Human Pathway Maps and Cellular Process Networks from the GeneGO database. The top-ranked pathways/processes (p < 0.05) included the endoplasmic reticulum (ER) stress response, neurophysiological process/visual perception, bile acid metabolism and negative FXR-dependent regulation of bile acid concentration, connective tissue degradation, and sodium transport (Figure 5E and Supplementary Data 1). Network visualization diagrams show upregulation of different branches of the ER stress pathway (Supplementary Figure 6), particularly downstream of the ATF-6 and IRE1 signaling pathways that activate pexophagy (sqstm1) and the protein folding response (GRP78, DnaJB9, Endoplasmin, hspa5). This indicates that peroxisomal dysfunction leads to ER stress and peroxisomal degradation. Visualization of the “neurophysiological process_visual perception” network highlighted downregulation of genes encoding especially light-sensitive receptors and alterations in their downstream pathways (Supplementary Figure 7). Interestingly, results showed a mild, non-significant decrease in all pex genes, suggesting that there is no transcriptional regulation of peroxisomal genes to compensate for Pex1 deficiency. The transcript levels of the pex1 gene itself were also only decreased by about 25% in the pex1^–/– larvae (as confirmed by qPCR; not shown), indicating that the mutant allele (17 bp deletion) does not undergo nonsense-mediated mRNA decay.

To assess potential neurological phenotypes, we analyzed the behavior of the larvae. First, we performed locomotion assays at different developmental stages and under normal light conditions. At 7 dpf, no significant differences were observed between the larvae of different genotypes (Figure 6A and Supplementary Figure 8A). However, starting at 11 dpf, pex1^–/– larvae showed a consistent decrease in movement, which became more evident at 13 and 15 dpf, indicating that their locomotion is progressively impaired in comparison to WT controls (Figure 6A and Supplementary Figure 8A). The difference in the locomotor behavior does not seem to be related to any obvious developmental defect, since pex1^–/– larvae show similar morphology and size compared to their WT siblings (Supplementary Figures 8B, C). To further evaluate locomotion-based behaviors, we performed startle response assays by delivering drastic light changes and tapping as stimuli. Similarly to the first assay, 7 dpf pex1^–/– and WT larvae showed no obvious differences in the expected behavioral response that is characterized by a transient period of hyperactivity as soon as there is a dark flash, followed by a progressive return to basal activity, which drops as soon as there is a light flash (Supplementary Figure 8D; Burgess and Granato, 2007). However, at later stages, pex1^–/– larvae displayed a shorter adaptation period to dark conditions than their WT sibling larvae, and the mutants did not show the sudden drop-in locomotor activity upon the light flash (Figure 6B). Interestingly, we did not observe a difference in the startle response upon tapping, which is an acoustic-vibrational stimulus generated by a built-in solenoid in the behavioral chamber and capable of inducing fast motor escape responses (Figure 6C and Supplementary Figure 8E). Such responses are mediated by Mauthner cells and associated neural pathways in the hindbrain (Liu and Fetcho, 1999).

The observations in the behavioral assays together with the transcriptomics data support the idea that visual perception pathways might be affected in the pex1^–/– larvae. Considering that ZSD patients develop retinal degeneration and glaucoma (Yergeau et al., 2023), we sought to assess whether similar hallmarks are present in our pex1^–/– zebrafish model. Eyes of 7-month-old pex1^–/– fish were cryosectioned and stained with DAPI. The inner nuclear layer (INL) did not show a significant difference in thickness between the pex1^–/– and WT samples (34.43 ± 1.82 μm and 34.49 ± 1.94 μm) (Figures 6D, E). However, the thickness of the outer nuclear layer (ONL) was significantly reduced in the pex1^–/– samples compared to the WT controls (14.01 ± 4.89 μm and 19.57 ± 6.23 μm, respectively) (Figures 6D, F). To further investigate changes in these layers, we performed immunohistochemistry on specific cell types. These analyses did not show statistically significant alterations in the overall number and structural organization (thickness) of Müller glial cells, spanning from the outer limiting membrane (OLM) to the ganglion cell layer (GCL) based on glutamine synthetase (GS) staining (Figures 6G–I). On the other hand, cone arrestin staining revealed a disorganized cone outer segment and a trend toward fewer cone cells compared to WT (Figures 6G, J). However, since the reduction in cone cell count was not statistically significant, the substantial thinning of the ONL suggests greater rod loss, leading to significant photoreceptor degeneration.