Ddhd1 knockout mouse as a model of locomotive and physiological abnormality in familial spastic paraplegia

Abstract We have previously reported a novel homozygous 4-bp deletion in DDHD1 as the responsible variant for spastic paraplegia type 28 (SPG28; OMIM#609340). The variant causes a frameshift, resulting in a functionally null allele in the patient. DDHD1 encodes phospholipase A1 (PLA1) catalyzing phosphatidylinositol to lysophosphatidylinositol (LPI). To clarify the pathogenic mechanism of SPG28, we established Ddhd1 knockout mice (Ddhd1[−/−]) carrying a 5-bp deletion in Ddhd1, resulting in a premature termination of translation at a position similar to that of the patient. We observed a significant decrease in foot–base angle (FBA) in aged Ddhd1(−/−) (24 months of age) and a significant decrease in LPI 20:4 (sn-2) in Ddhd1(−/−) cerebra (26 months of age). These changes in FBA were not observed in 14 months of age. We also observed significant changes of expression levels of 22 genes in the Ddhd1(−/−) cerebra (26 months of age). Gene Ontology (GO) terms relating to the nervous system and cell–cell communications were significantly enriched. We conclude that the reduced signaling of LPI 20:4 (sn-2) by PLA1 dysfunction is responsible for the locomotive abnormality in SPG28, further suggesting that the reduction of downstream signaling such as GPR55 which is agonized by LPI is involved in the pathogenesis of SPG28.


Introduction
Spastic paraplegias (SPGs) are neurological disorders characterized by spasticity and gait disturbance. Abnormalities of the pyramidal tract are known to be a hallmark of this disorder. More than 60 types of SPGs have been reported to be genetically distinct [1]. SPG type 28 (SPG28) is an autosomal recessive SPG caused by mutations in the gene encoding DDHD domain-containing protein 1 (DDHD1) also known as phospholipase A 1 (PLA 1 ) [1][2][3][4][5]. We have previously identified a novel homozygous 4-bp deletion (c.914 917delGTAA, p.Ser 305 Ilefs*2) in exon 2 of the DDHD1 gene as the variant responsible for SPG28 (OMIM#609324) [6]. Phospholipase A is known to catalyze phosphatidylinositol (PI) and phosphatidic acid (PA) to lysophosphatidylinositol (LPI) and lysophosphatidic acid (LPA), respectively. There are two kinds of phospholipase A; PLA 1 such as DDHD1, hydrolyzes the sn-1 ester bond of PI, and phospholipase A 2 (PLA 2 ) hydrolyzes the sn-2 ester bond of PI [7]. The intracellular PLA 1 protein family is characterized by the presence of a short lipase active-site sequence and a C-terminal DDHD domain. DDHD1 is known to be highly expressed in the brain and testes and its dysfunction causes neurodegeneration with brain iron accumulation (NBIA) as well as SPG phenotypes [8][9]. Animal models are useful in clarification the pathogenic mechanism and potential therapy. There are two previous studies of Ddhd1 knockout (KO) mice [8,10]. One has reported impaired movement of sperm and abnormal mitochondrial morphology with no description of abnormal mobility [8]. The other study has reported a significant increase in PI 18:1/20:4 and a significant decrease in LPI 20:4 in cerebra of their Ddhd1 KO mice. The authors, however have not described any abnormal mobility in 6-month-old Ddhd1 KO mice [10]. SPG onset in humans varies at ages from 0 to 70 years and deteriorates in aged patients, suggesting that the mice examined in these previous studies might have been too young for the authors to examine SPG phenotypes [2]. To clarify the pathogenic mechanism of SPG28, we established Ddhd1(−/−) and perform behavioral analyses (14 and 24 months of age). We also performed RNA sequencing and lipidome analysis on sufficiently aged mice (26 months of age).

Establishment of Ddhd1 KO mice
We generated Ddhd1 KO mice using the CRISPR/Cas9 system according to Yang et al. (2014) [11] ( Figure 1A). We identified four kinds of small indels that are expected to cause frameshift and premature termination at Ddhd1. We selected a mouse harboring 5-bp deletion in exon 2 since the resulting amino acid sequence is the closest to the variant found in SPG28 patients among some kinds of indels created by CRISPR/Cas9 ( Figure 1B,C) [6]. To remove off-target sites from the mice harboring 5-bp deletion, we performed one-generation backcross with C57BL/6J. By crossing the F 1 mice heterozygous for the 5-bp deletion, we established the Ddhd1(−/−) strain as the Ddhd1 KO mice carrying the premature termination at a very similar position with the variant we identified in the original patient. We also confirmed the absence of DDHD1 protein signal in Ddhd1(−/−) mice by Western blotting ( Figure  1D). Full-length blots are presented in Supplementary Figure S1. In the following analyses, we use heterozygous KO mice (Dhdh1[+/−]) as a control for homozygous KO mice (Ddhd1[−/−]).

RNA sequencing and gene ontology enrichment analysis
We sequenced total RNA extracted from the cerebra of two Ddhd1(+/−) and two Ddhd1(−/−). We obtained an average of 30609158 reads per mouse used in the present study.  Table S4). Using gene ontology enrichment analysis (GEA), we identified two annotation clusters which include significant gene ontology (GO) terms (Benjamini probabilities < 0.05) as shown in Table 1. We identified significantly enriched terms related to the nervous system as GO:0007268 (synaptic transmission), GO:0051966 (regulation of synaptic transmission, glutamatergic), GO:0019226 (transmission of nerve impulse), and GO:0007270 (nerve-nerve synaptic transmission). GO terms related to cell-cell communication, such as GO:0007267 (cell-cell signaling), were also enriched as some of significant GO terms (Table  1). We have previously reported that the decrease in DDHD1 mRNA expression level in the peripheral blood of SPG28 patients compared with unrelated control due to the nonsense-mediated decay [6]. Although Ddhd1 was not selected as a significant DEG in the current RNA-seq analysis, we confirmed the decrease in the Ddhd1 transcription level by approximately 55% in the Ddhd1(−/−) mice cerebra (Supplementary Table S5). Among the DEGs, we examined the expression of Rtn4r by real-time quantitative PCR (RT-qPCR) analysis. We observed significant increase in Rtn4r and Adra2a expression in Ddhd1(−/−) at the 12 months of age consistent with the result of RNA sequencing (Supplementary Figure S6).

Immunohistochemistry
We examined neurofilament (NF) protein expression as an axonal marker in spinal of mice at the 6 months of age. While there was no apparent morphological difference between Ddhd1(−/−) and Ddhd1(+/−), we observed remarkable decrease in staining of NF in the pyramidal tract of Ddhd1(−/−) compared with Ddhd1(+/−) ( Figure  5A-C). This result indicates axonal decrease in the pyramidal tract in Ddhd1(−/−) at 6 months of age.

Discussion
In the present study, we observed that gait disturbance in SPG patients was partially replicated in our Ddhd1 KO mice, Ddhd1(−/−) at 24 months of age. Although there are two previous reports of Ddhd1 KO mice [8,10], neither has reported locomotive abnormality. This discrepancy can be attributed to the difference in ages of mice examined. In one of the previous studies, locomotive phenotypes of Ddhd1 KO mice were examined up to 6 months of age [10], while we examined FBA in 24 months of age, which are the oldest Ddhd1 KO mice ever examined for their locomotion. Consistently we failed to observe the FBA phenotype in younger mice at the 14 months of age. This suggests that SPG phenotypes takes long time to manifest such as at least more than 14 months for the FBA phenotypes in Ddhd1 KO mice. Fourteen months of age in mice is equivalent to forties in humans [15]. It is consistent with the observation that the age of SPG onset in human is known to be variable and sometimes to be seventies [2]. Therefore, even longer observation (>24 months) can be helpful to detect the SPG phenotypes in milder SPG mice. It is quite possible that abnormal locomotion could have been observed in other Ddhd1 KO mice previously reported through the long-term observation. In addition, examination of the expression levels of genes involved in the LPI metabolism  pathway in different ages may be helpful to clarify the mechanism of the deterioration of the phenotype along with aging.
The FBA test was originally established as an approach to assess muscle function in mice that had experienced femoral nerve damage [16]. Even though FBA is one of the simplest methods to quantitate locomotive activity without requiring special devices, it has not been applied to the examination of Ddhd1 KO mice. Two strains of SPG model mice (subtypes SPG15 and SPG31) have also been reported to show significant decreases in their FBA [12,13], suggesting that FBA is an easy and reliable method in quantitative phenotyping of locomotion. From these studies, SPG15 and SPG31 model mice have observed symptoms by measuring FBA at 12 and 4 months, respectively. In contrast, we showed SPG28 takes much longer time to manifest the symptoms. As ages of SPG onset is known to be variable in humans, it is reasonable that ages of onset in SPG model mice are also variable among the types of SPGs.
Although gait disturbance is the common symptom in spinocerebellar ataxia (SCA) as well as in SPGs [23], their signs can be clearly distinguished by neurological examination. DDHD1 proteins are highly expressed in the cerebella as well as in the cerebra (Supplementary Figure S7). However, our previous neurological examination has confirmed no signs of cerebellar ataxia in the SPG28 patient in the original pedigree [6], indicating there is no involvement of cerebellum in the gait disturbance we observed in Ddhd1 KO mice. Interestingly, the amount of LPI 20:4 (sn-2) decreased in the cerebra, but did not significantly change in the cerebella (P=0.09) (Supplementary Figure S4). While hindpaws of SCA model mice have been reported to be slipped off of the beam [24], we did not observe such slip-off in our Ddhd1 KO mice during beam-walk test (data not shown). These data suggest that symptoms of gait disturbance is caused by a decrease in LPI in the cerebra, not in the cerebella.
In the present study we observed significant differences were observed in multiple PIs such as PI 16 Some GO terms significantly enriched in DEGs in Ddhd1(−/−) mouse. Two annotation clusters included significant GO terms (Benjamini probability < 0.05). GO terms related to nervous system or cell-cell communications are marked with asterisks (*).
PI 18:1/20:4 [10]. This discrepancy can also be explained by the difference of examined mice (26 months in our study and 1.5-3 months in [10]). Although the previous study has reported the decrease in LPI 20:4 in Ddhd1 KO mice [10], the researchers examined only total amount of LPIs without distinguishing their components, such as LPI 20:4 (sn-1) or LPI 20:4 (sn-2). In the current study, we observed the specific decrease in LPI 20:4 (sn-2) rather than that of LPI 20:4 (sn-1) in Ddhd1(−/−) by distinguishing isoforms of LPIs using SFC. Our result suggests that the decrease in total LPI 20:4 observed in previous study is due to a decrease in LPI 20:4 (sn-2), not LPI 20:4 (sn-1). LPIs are known to be endogenous agonists of GPR55, which might be involved in the regulation of axon growth and synaptic formation [17][18][19][20][21][22]. This observation suggests that the reduction in LPI 20:4 (sn-2) suppresses the signaling of the GPR55 and triggers abnormal axonal extension and incomplete neural circuits, eventually resulting in abnormal locomotion. To further examine the pathogenic mechanism, it is necessary to perform proteome analyses and functional analyses of neural cell formation in Ddhd1(−/−). PA consumption and LPA production by the PA-PLA 1 activity have been suggested to play an important role for mitochondrial fission [10]. However, our lipidome analyses observed no significant changes of PAs nor LPAs amount in Ddhd1(−/−) (Supplementary Table S3), suggesting that mitochondrial phenotypes may be associated with specific types of PIs and LPIs rather than PAs and LPAs.
In our current analysis, only three genes show significantly increased expression in cerebral tissues of Ddhd1(−/−), although we did not directly examine their protein expression levels. One of these, Rtn4r, is notable since it is known to mediate axonal growth inhibition and is involved in the regulation and plasticity of the adult central nervous system [25]. Since Rtn4r is known to activate RhoA, which is a downstream molecule of the GPR55 signaling, the increased expression of Rtn4r is possibly a compensatory feedback to recover the abnormal deactivation of GPR55 caused by Ddhd1 dysfunction [26]. Interestingly, three genes encoding G protein-coupled receptors (GPCRs) were listed in DEGs, Adra2a, Drd1a and Drd2 (Table 1 and Supplementary Table S4). Although these GPCRs do not interact with LPIs, the decrease in GPR55 signaling may disturb intracellular GTP metabolism, resulting in decreased expression of these GPCR genes.
We observed axonal degradation in the pyramidal tract in Ddhd1(−/−) at 6 months of age, which is much younger than the age we observed the weakness in hindlimbs by FBA analysis (24 months) ( Figure 5). Therefore, axonal degradation in the pyramidal tract is likely to be one of the preclinical neurological changes preceding the onset of SPG phenotype. These results suggest that our Ddhd1 KO mice at the age of 6 months are preclinical stage of SPG28 only showing a decrease in the number of axons.
DDHD2, another gene encoding PLA 1 , is also known to be responsible for SPG54 [27,28]. In addition, two other genes involved in the LPI metabolism pathway, NTE and CYP2U1, are also known to be responsible for SPG39 and SPG56, respectively [1,3]. Therefore, the disruption of the normal metabolism of LPIs is likely to be the common pathogenic mechanism shared by multiple SPGs, such as SPG28, SPG39, SPG54, and SPG56, suggesting the potential common drug for multiple SPGs.
SPG model mice reported here provide two major potential applications. One is the fine phenotyping of the disease progression. Physiological changes observed in the model mice provide seeds of biomarkers at the subclinical stages of SPG. The second is the platform for the development of therapeutic drugs. Candidate molecules can be deduced from the molecules connected to the DEGs identified in Ddhd1(−/−). Model mice can be examined at arbitrary ages so that the disease progression, the drug effects and the kinetics of candidate biomarkers can be analyzed through chronological snapshots.

Materials and methods
Animals B6C3F1 and ICR mice were obtained from Kyudo company (Saga, Japan), and C57BL/6J mice were obtained from Charles River Laboratories Japan (Yokohama, Japan). They were fed a standard pellet diet (CLEA Japan, Inc., Tokyo, Japan) and filtered water. The animals were kept under condition of a 12:12-h light:dark cycle. All animal experiments were performed in Medical Institute of Bioregulation, Kyushu University.

CRISPR construct generation
px330-U6-Chimeric BB-CBh-hSpCas9 plasmid vectors harboring the cDNA sequence encoding Streptococcus pyogenes Cas9 (hSp-Cas9) and AmpR was purchased from Addgene (catalog 42230) [29]. Guide sequence oligonucleotides Mouse CRISPR target S and Mouse CRISPR target AS are shown in Supplementary Table S6. Oligonucleotides were annealed (95 • C for 10 min followed by cooling down at room temperature for 30 min) and cloned in the px330 vector, which was digested with BbsI. Constructs were introduced into competent DH5α cells. The colonies harboring relevant constructs were inoculated into 5 ml Luria-Bertani medium and cultured overnight. Plasmid DNA was extracted from the culture using a Plasmid Maxi kit (#12165, Qiagen, Hilden, Germany).

Microinjection
Pregnant mare serum gonadotropin and human chorionic gonadotropin were injected into B6C3H female mice at 48-h intervals. Injected females were then housed with C57BL/6 male mice. Fertilized one-celled stage embryos were collected from oviducts. We performed zygote injection according to Yang et al. (2014) [12]. We held a zygote using a holding pipette, inserted the injection pipette into the zygote without breaking the oolemma, and advanced the pipette until it almost reached the opposite side of the zygote's cortex. Approximately 5 pg of px330 was injected into each zygote. Injected zygotes were cultured in KSOM medium (#MR-121-D, MERCK, Darmstadt, Germany) at 37 • C in a 5% CO 2 incubator until the two-celled stage and were then transferred into pseudo-pregnant ICR mice to obtain the initial generation of genome-edited mice, F 0 . The protocols were approved by the Institutional Animal Care and Use Committee of Kyushu University.

Establishment of the mouse strain
The F 0 generation mice determined to be positive by the surveyor nuclease assay were bred with wildtype C57BL/6J to produce the F 1 mice, which were identified by PCR and sequencing. Male infertility has been suggested in

DNA preparation
We collected tails from 2 weeks of age mice for genotyping. The tails were immersed in 50 μM NaOH at 95 • C for 10 min and centrifuged at 12000 rpm for 15 min. We performed PCR using the supernatant fluid as a PCR template.

Sanger sequencing
Exon 2 of Ddhd1 was sequenced by direct sequencing using the primers Mouse Ddhd1 forward and Mouse Ddhd1 reverse (Supplementary Table S7

FBA analysis
We allowed mice to walk on an elevated horizontal plastic beam (50 × 5 cm) with 20 cm of the length and shot their videos from behind by iPhone (Apple, Cupertino, CA, U.S.A.). Once the mouse walked to the middle of plastic beam, it was brought back to the front with its tail until more than 20 steps were recorded. Angles at the toe-off positions of hind paws of was measured as FBA for the first 20 steps using single video frames at 14 months of age (n=3) and 24 months of age (n=2) awere examined. The data between Ddhd1(+/−) and Ddhd1(−/−) were analyzed using a two-tailed Student's t test and the data between individuals were analyzed using a Tukey-Kramer test. Tukey-Kramer test was performed in R Studio (version 3.5.3).

Animals used in tissue sampling
At 26 months of age, Ddhd1(+/−) and Ddhd1(−/−) (n=2 per group) were killed with cervical spine fracture dislocation without anesthesia. Cerebra were extracted from each mouse and immediately snap-frozen in liquid nitrogen for subsequent lipid and protein isolation. The cerebra samples for RNA-seq were immediately immersed in RNAlater Solutions (#AM7021, Invitrogen, Carlsbad, CA, U.S.A.). We utilized all available animals (two mice each for Ddhd1(+/−) and Ddhd1(−/−)) after the long-term rearing of 26 months. Details of the mice used in FBA tests, lipidome analyses and RNA sequencing were described in Supplementary Table S8.

Western blotting
Since the cerebrum has been reported to show the highest Ddhd1 expression [8], we extracted total protein from cere-

RNA sequencing
Total RNAs were extracted from the cerebra of both 26-month-old Ddhd1(+/−) and Ddhd1(−/−) mice (n=2 per group) using the RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol. cDNA libraries were constructed from the purified RNA using the SMARTer Universal Low Input RNA Kit for Sequencing (#634938, Clontech, Mountain View, CA) and a Low Input Library Prep Kit (v 2) (#634947, Clontech, Mountain View, CA). Single-read 50-bp sequencing was performed on the Illumina HiSeq 2500. To generate a heatmap, we used the 'heatmap.2' function of the gplots package in R software. This analysis was performed in R Studio (version 3.5.3).

Detection of DEGs and GEA
Raw reads obtained from NGS were mapped to the reference genome (Mouse GRCm38/mm19) using TopHat. DEGs were identified using Cuffdiff to test the significance of the differential expression of genes based on fragments per kilobase of exon per million fragments mapped (FPKM). Finally, The DEGs were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.7, https://david.ncifcrf.gov/home.jsp) to analyze the enrichment of GO terms. GO terms showing the Benjamini probabilities < 0.05 were considered as significant difference.

Real-time quantity PCR
The purified RNA extracted from the cerebra of both 12-month-old Ddhd1(+/−) and Ddhd1(−/−) mice (n=2 per group) was converted in each instance into cDNA using ReverTra Ace qPCR RT Master Mix with gDNA Remover (TOYOBO, Osaka, Japan). RT-qPCR were then performed with Power SYBR Green Master Mix (#4367659, Applied Biosystems, Waltham, MA) on the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, Waltham, MA). Expression of Gapdh was also examined as an internal standard of mRNA expression. Primers are shown in Supplementary Table S7.