Zebrafish of the AB strain used for the experiments were maintained in the fish facilities of the Department of Zoology, Genetics and Physical Anthropology of the University of Santiago de Compostela. They were maintained at 28°C with a photoperiod of 14 h of light and 10 h of darkness according to previously described protocols (Westerfield, 2000; Aleström et al., 2019). All experiments involving animals followed the guidelines of the European Community and Spanish Government on animal care and experimentation (Directive 2012-63-UE and RD 53/2013) and were approved by the bioethics committee of the University de Santiago de Compostela and the Xunta de Galicia government.

Zebrafish npc1 shows 70% sequence identity with its human ortholog. Zebrafish npc1 (Ensembl ID: ENSDARG00000017180) sequence is composed of 25 exons and codifies for 1,276 amino acids. We created the mutation at the beginning of the exon 22 because this is part of the cysteine-rich luminal loop where the most common sites for mutations in NPC1 related to juvenile and severe onset of the disease are located (Millat et al., 2001; Ribeiro et al., 2001; Park et al., 2003; Gelsthorpe et al., 2008).

Two guide RNAs (gRNAs) placed at exon 22 (which shares 78% sequence identity with the human exon 22) of the npc1 zebrafish gene were designed with CRISPRscan¹:

gRNA1: taatacgactcactataGGAAGACGTAAAACACACTGgttttagagctagaa.

gRNA2: taatacgactcactataGGGCTTTGAGCTCTGGTCGGgttttagagctagaa.

Each gRNA primer was amplified with the universal primer 5′-AAAAGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACGGACTAGCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC-3′, using iProof^™ High-Fidelity DNA Polymerase (Bio-Rad; Hercules, CA, United States) (protocol available upon request). Subsequently, gRNAs were in vitro transcribed with MAXIscript^™ T7 Transcription Kit (Thermofisher Scientific; Waltham, MA, United States) and gRNA concentration was measured through NanoDrop^® 2000 (Thermo Fisher Scientific) spectrophotometer.

Between 200 and 300 zebrafish embryos at one cell stage were injected with a mix of 15–40 ng/μL of each gRNA and 1.2 μg/μL of TrueCut^™ Cas9 Protein v2 (Thermofisher Scientific). Each embryo was injected with 5 nL of this mixture. The experiment was repeated three times.

To verify efficiency of the CRISPR/Cas9 protocol, genomic DNA from 10 injected 48 h post-fertilization (hpf) embryos was extracted with Chelex 100 Resin (Bio-Rad) and PCR amplified with AmpliTaq Gold^™ DNA Polymerase (Thermo Fisher Scientific) (protocol available upon request) using the following primers npc1-F: 5′-GGACTTCCCAGTGACATAATGG-3′ and npc1-R: 5′-CCTGGTGCTGATGGAGAAAG-3. The presence of the heteroduplex product of the CRISPR/Cas9 editing was checked in a 5% acrylamide gel.

F0 mutant embryos were raised to adulthood and out-crossed with wt individuals. The resulting F1 generation was also raised to adulthood and checked for CRISPR/Cas9 mutations. Heterozygous npc1 mutant individuals (around 50% of the fish) were distinguished from wt through genotyping DNA obtained from the caudal fin. Amplified DNA from F1 mutants was used for cloning with the StrataClone PCR Cloning Kit (Agilent; Santa Clara, CA, United States) and then transformed with StrataClone SoloPack (Agilent) in LB-ampicillin plates. Thereafter, 10 colonies of each npc1 heterozygous mutant were amplified and double checked by 1% agarose gel electrophoresis and Sanger sequencing with 3730×l DNA Analyzer (Thermo Fisher Scientific). Finally, npc1 heterozygous zebrafish (npc1^+/−) of the F1 carrying the same mutation were in-crossed to obtain F2 homozygous npc1 (npc1^−/−), which were used for subsequent experiments. All analyses were performed without prior knowledge of the genotype of the individuals, with the exception of RNAseq and lipid analysis in which we needed to know the genotype in order to create the pools of larvae for shipping to the service providers.

1 to 3 weeks post-fertilization (wpf) F2 fish from npc1^Δ56 and npc1^Δ7 lines carrying wt, heterozygous (npc1^+/−) and homozygous (npc1^−/−) mutant alleles were used for phenotypic characterization. For motor behavior analysis, body length analysis and histopathological analyses, we have used the same larvae (for histopathological analyses only 2 and 3 wpf larvae) resulting in:

For this purpose, larvae were anesthetized with 0.002% tricaine methanesulfonate (MS-222, Sigma-Aldrich; Saint Louis, MO, United States). Photographs were taken every week with a Nikon Ds-Ri1 camera attached to an inverted fluorescence microscope (AZ100 Multizoom Nikon; Tokyo, Japan). Images were analyzed with ImageJ software (National Institutes of Health; Bethesda, MD, United States) (Abràmoff et al., 2004) and total body length was measured. All fish were genotyped at the end of the experiments.

For immunohistochemical analysis, 3 μm sections of F2 fish of 3 wpf, 5 fish per genotype (wt, npc1^+/−, npc1^−/−) of npc1^Δ56 and npc1^Δ7 lines were dewaxed and rehydrated at room temperature in a humid chamber. Slides were washed with 10 mM phosphate-buffered saline and 0.5% Tween 20 (pH 7.4) in three successive 5 min immersions. Endogenous peroxidase was quenched with Bloxall Blocking Solution (Vector Laboratories, Inc.; Burlingame, CA, United States) and the sections were incubated with a primary rabbit monoclonal recombinant anti-Niemann Pick C1 antibody [EPR5209] (ab134113 Abcam; Cambridge, United Kingdom) generated against the C-terminal region of the protein. Immunohistochemical staining was carried out with ImmPRESS HRP Horse Anti-Rabbit IgG Polymer Kit (Vector Laboratories, Inc.), together with the Vector VIP Substrate Kit (Vector Laboratories, Inc.) as chromogen. Negative controls were included in the experiment by using the antibody dilution solution without primary antibody. Sections were counterstained with hematoxylin.

For survival analysis, embryos resulting from the cross of heterozygous adults were used. Embryos and larvae were maintained for the desired time, by the end of which surviving fish were euthanized and genotyped. The number of individuals in each group at the start of the experiments were: 1 wpf: 179 individuals: 86 for npc1^+/−Δ56 intercross and 93 for npc1^+/−Δ7 intercross; 2 wpf: 257 individuals: 163 for npc1^+/−Δ56 intercross and 94 for npc1^+/−Δ7 intercross, 3wpf: 188 individuals: 126 for npc1^+/−Δ56 intercross and 62 for npc1^+/−Δ7 intercross; 4 wpf: 86 individuals: 44 for npc1^+/−Δ56 intercross and 42 for npc1^+/−Δ7 intercross; and 5 wpf: 40 individuals: 19 for npc1^+/−Δ56 intercross and 21 for npc1^+/−Δ7 intercross. Number of surviving fish for each genotype were annotated and represented using a clustered bar chart of percentages.

Motor performance of F2 zebrafish larvae (npc1^Δ56 and npc1^Δ7) at 1, 2, and 3 wpf was quantified with a Zebralab system with a Zebrabox (Viewpoint; Civrieux, France). The quantification software measures the number of pixels moved for each fish in a certain period. For this, 1 wpf larvae were introduced in 96-well plates, while 2 and 3 wpf larvae were placed in 24-well plates. Total larvae movement was measured for 1 hour, alternating 10 min periods of light vs dark conditions.

Statistical analysis and graphs of the morphology and locomotion of the three groups (wt, npc1^+/− and npc1^−/−) were generated with GraphPad Prism version 7 (GraphPad; San Diego, CA, United States), using a non-parametric Kruskal-Wallis test with Dunn’s multiple comparisons test after performing the D’Agostino–Pearson normality test. Statistical significance was established at a value of p < 0.05. Mean with Standard Deviation (SD) are represented in the graphs.

Sphingomyelin bodipy (D3522, Thermo Fisher Scientific) and Topfluor Cholesterol (810255P, Avanti; Webster, United States) stock solutions were produced by dissolving 1 mg in 2 mL chloroform, and 20 μL of stock solution was transferred to eppendorf tubes to let evaporate. The powder was then dissolved in DMSO at a concentration of 10 mg/L. In vivo staining of 1 to 2 wpf npc1^Δ56 larvae was carried out for 24 h at a final concentration of 0.01 mg/L. After staining, larvae were rinsed 3 times with fish water and anesthetized with 0.002% tricaine methanesulfonate.

Confocal photomicrographs were taken with a Leica TCS SPE confocal laser microscope (Leica Microsystems; Wetzlar, Germany). After capturing the images, larvae were used to verify the genotype.

The images were analyzed with ImageJ software (National Institutes of Health; Bethesda, MD, United States). Pixels of the embryo were selected using the tracing tool and the area, the mean and the integrated density of each fish were measured. These measurements were combined in the formula of corrected total cell fluorescence (CTCF): CTCF = integrated density − (area of selected cell × mean fluorescence of background readings). Statistical analysis and graphs were generated with GraphPad Prism version 7, using a Mann–Whitney test after performing the D’Agostino–Pearson normality test. Statistical significance was established at a value of p < 0.05. Mean with SD are represented in the graphs.

For RNAseq analyses, 6 samples were used (3 wt and 3 npc1^{Δ56/ Δ56}) each of them containing a pool of 15–20 2 wpf F2 larvae, resulting from F1 incrosses. After sampling the end of the tail for genotyping, larvae were immersed in RNAlater (Sigma-Aldrich) and conserved at −20°C. RNA extraction was performed with miRNeasy Micro Kit (Qiagen; Hilden, Germany). RNA was quantified using NanoDrop^® 2000 (Thermo Fisher Scientific), and a total of ~1,000 ng (40 ng/μL) was sent for sequencing. RNA integrity was evaluated in a 2100 Bioanalyzer (Agilent Technologies), and all samples displayed RNA integrity number (RIN) values >9 so appropriate for library construction and sequencing. Libraries were sequenced on an Illumina Nova-Seq 150 bp PE run by Novogene (Cambridge, United Kingdon). The quality of the sequencing output was assessed using FastQC v.0.11.7.²

Removal of low quality reads, contaminating sequence, low complexity reads (repetitive DNA sequences) and short reads was performed on read pairs using Fastp v.0.20.1 (Chen et al., 2018). Illumina specific adaptors as well as leading and trailing bases with a Phred score < 20 were eliminated; only reads where both pairs were longer than 35 bp post-filtering were retained.

Alignment of filtered reads against the latest version of the zebrafish genome (GRCz11)³ and quantification of transcript abundance was performed using kallisto v.0.46.1 (Bray et al., 2016). According to the authors’ description, Kallisto “pseudo aligns” reads to the transcriptome, producing a list of transcripts that are compatible with each read. This software avoids alignment of individual bases against the genome, thus allowing accounting for multi-mapping reads. Aligned reads were assigned to genes based on the latest annotation of the zebrafish genome.

Gene count data were used to calculate gene expression and estimate differential expression (DE) between wt and mutants, using the Bioconductor package DESeq2 v.1.28.1 (Love et al., 2014) in Rstudio v.4.1.2 (R Core Team, 2017). Briefly, size factors were calculated for each sample using the ‘median of ratios’ method and count data were normalized to account for differences in library depth. Next, gene-wise dispersion estimates of gene counts were fitted to the mean intensity using a parametric model and reduced toward the expected dispersion values. Differential gene expression was evaluated using a negative binomial model fitted for each gene, and the significance of the coefficients was assessed using the Wald test. The Benjamini-Hochberg false discovery rate (FDR) correction for multiple tests was applied, and transcripts with FDR < 0.05 were considered differentially expressed genes (DEGs). Hierarchical clustering and principal component analyses (PCA) were performed to assess sample clustering and identify potential outliers over the general gene expression background. Gene Ontology (GO) enrichment analysis of the DEGs was carried out with ShinyGO v0.741 (Ge et al., 2020), using the zebrafish transcriptome as background.

Tissue for lipid profiling was obtained from F2 larvae of 2 wpf of npc1^Δ56. All individuals were sectioned into two parts, most of the body and the end of the tail. The body was frozen at −80°C and the end of the tails were used to determinate the genotype. 5 pools of 20–25 wt and 5 pools of 20–25 npc1^{Δ56/ Δ56} larvae were employed for the analysis.

The 10 samples were sent to Lipotype GmbH (Dresden, Germany) for lipid profiling. Concentrations of phosphatidate, phosphatidylcholine, phosphatidylethanolamine, phosphatidylglycerol, phosphatidylinositol, phosphatidylserine, diacylglycerol, triacylglycerol, sphingomyelin and cholesteryl ester were measured by a mass spectrometer and a final bioinformatic data report was generated by the same company (Surma et al., 2015).

In summary, the amount of a lipid class was calculated by summing the pmol values of the individual lipids belonging to each class. Class amount is then normalized to total lipid content. For normalization, internal standards (one per class) were used. The intensities of the endogenous lipids are normalized to the intensities of the respective standard and based on the used amount of the standard the absolute value of the measured lipid can be calculated. Values are mean of the biological replicates. Statistical significance was calculated by ANOVA after performing the D’Agostino–Pearson normality test. Statistical significance was established at a value of p < 0.05. Mean with SD are represented in the graphs.

We verified by acrylamide gel electrophoresis the presence of an heteroduplex product which indicates the efficiency of the CRISPR/Cas9 technique. Two different mutant lines of npc1 were identified:

Then, the npc1^Δ56 line was genotyped by agarose gel electrophoresis, while the npc1^Δ7 line was genotyped by fragment analysis on an automated fluorescent sequencer (3730xl DNA Analyzer). Both mutations were placed at the C-terminal part of npc1, at the end of the cysteine rich luminal loop between the transmembrane domains 8 and 9 and were predicted to cause a premature STOP codon (Figure 1).

The preliminary characterization of the homozygous mutants (npc1^{Δ56/ Δ56} and npc1^{Δ7/ Δ7}) revealed similarities in histopathology, survival and locomotion in both lines and therefore these analyses are reported together for both lines by referring to them as npc1^−/− (see below). The RNAseq and the lipid staining and lipid profiling were only performed in the npc1^Δ56 line, due to its greater ease of genotyping.

Npc1^−/− and wt individuals showed no apparent macroscopical and histological differences until 2 wpf, albeit only 80% of npc1^−/− larvae reached this age. The mortality rate was higher at 3 wpf with a survival of 60%. At 4 wpf survival rate of npc1^−/− larvae was 12% and none of them survived up to 5 wpf (Figure 2). From 2 wpf, npc1^−/− showed a significantly shorter body length (value of p = 0.0343; mean length: wt: 4,155 μm ± (SD): 717; npc1^+/−: 3,998 μm ± 595.3; npc1^−/−: 3,633 μm ± 640.8). From 3 wpf, the smaller size became much more evident (value of p = 0.0037; mean length: wt: 5,523 μm ± 889.8; npc1^+/−: 5,358 μm ± 1,096; npc1^−/−: 4,376 μm ± 696.5) (Figure 3).

A histopathological analysis in 2–3 wpf npc1^−/− larvae revealed well-defined and optically empty intracellular spherical vacuoles, which were compatible with the accumulation of lipids in the liver (Figure 4), intestine, renal tubules, and occasionally in the brain (Figure 5). Severity of these lesions varied between different npc1^−/− individuals, however the liver was the location of greatest damage in all studied fish. Histopathological grades of vacuolation were assigned as mild, moderate and severe based on an increasing extent of the lesion and the complexity of the changes in the analyzed tissues. The lesions that were found in npc1^−/− larvae were considered moderate to severe degree.

We then aimed to characterize the presence of NPC1 protein in the internal organs by immunohistochemical studies. Our results showed a complete lack of immunostaining in all organs of npc1^−/− larvae, whereas wt individuals showed immunostaining in the central nervous system, liver, intestine, and renal tubules (Figure 6).

Finally, both cholesterol and sphingomyelin staining resulted in the similar stainings with accumulations of those lipids mainly in the liver of npc1 ^−/− larvae but also in the intestine (Figure 7).

At 1 wpf, no differences were found in locomotion between wt, npc1^+/− and npc1^−/− larvae, as measured through comparison of number of pixels moved (not shown). However, at 2 (not shown) and 3 wpf, there was a statistically significant reduction in locomotion of npc1^−/− larvae in comparison with wild type and npc1^+/− (value of p = 0.0002; mean pixels moved: wt: 133,222 ± 142,726; npc1^+/−: 150,380 ± 173,416; npc1^−/−: 68,921 ± 64,309) (Figure 8).

RNAseq analysis showed 284 differentially expressed genes (DEGs) between npc1^Δ56/Δ56 and wt 2 wpf larvae, of which 99 show a higher expression in npc1^Δ56/Δ56 animals and 185 have a reduced expression in npc1^Δ56/Δ56 in comparison with wt (Supplementary Table 1). Samples of mutant vs control individuals did not clearly cluster apart in the principal component analysis (PCA), suggesting the presence of individual variability (Figure 9). However, since the npc1 gene is significantly downregulated in all the npc1^Δ56/Δ56 samples our results can be considered valid to study differential gene expression and perform enrichment analyses.

Gene ontology (GO) analysis revealed 357 enriched biological processes, of which 32 involved the npc1 gene, like ‘lipid homeostasis’ and ‘response to hypoxia.’ We also found 111 enriched molecular functions such as ‘transporter activity,’ in some of them npc1 was involved. Other important molecular functions were ‘haptoglobin binding’ and ‘heme binding’ (Supplementary Table 2).

Out of the 284 DEGs, 14 are involved in transmembrane transporter activity (Supplementary Table 1). Additionally, we identified 17 DEGs related to lipid transport, modification and metabolic processes; 31 DEGs related to processes of central nervous system such as neural plasticity, axonogenesis, synaptogenesis, myelination and neurotransmitter activity (Figure 10); and 12 DEGs related to muscle contraction and development. We also found 12 DEGs related to cytoskeleton activity (Supplementary Table 1).

Among all DEGs, 20 genes related to ‘collagen extracellular matrix structural constituent’ were down-regulated in npc1^Δ56/Δ56, Furthermore, 17 DEGs related to hematopoiesis and heme oxygenase activity and 9 genes related with angiogenesis were deregulated in npc1^{Δ56/ Δ56} (Supplementary Table 1).

Five of the 10 lipid classes of lipids measured showed statistically significant differences in concentrations between npc1^{Δ56/ Δ56} and wt (Figure 11).

Principal component analysis showed a tendency for sample clustering and separation but with variability among samples of the same group (Figure 12).

Cholesteryl esters showed an overall statistically significant reduction in npc1^{Δ56/ Δ56} larvae (value of p = 0.02) compared to wt. 5 cholesteryl species were reduced in npc1^{Δ56/ Δ56} while another 10 species detected in wt larvae could not be found in npc1^{Δ56/ Δ56} (Figure 12).

Sphingomyelins showed a statistically significant increase in npc1^{Δ56/ Δ56} samples (value of p = 0.002). 11 sphingomyelin species were found to be more abundant in npc1^{Δ56/ Δ56}larvae than in wt and another 5 species detected in npc1^{Δ56/ Δ56} could not be found in wt (Figure 13).

Some of the phospholipids analyzed showed statistically significant differences in concentration between wt and npc1^{Δ56/ Δ56} larvae pools. Phosphatidate had a higher abundance in wt than in npc1^{Δ56/ Δ56}(value of p = 0.036). Phosphatidylglycerol and phosphatidylinositol, were reduced in wt compared with npc1^{Δ56/ Δ56} (value of p = 0.0055 and value of p = 0.038, respectively). Concentrations of other lipids (diacylglycerol, triacylglycerol phosphatidylcholine, phosphatidylethanolamine, and phosphatidylserine) did not show significant differences between wt and npc1^{Δ56/ Δ56} larvae.