Abstract Background Many psychiatric disorders share a working memory (WM) impairment phenotype, yet the genetic causes remain unclear. Here, we generated genetic profiles of WM deficits using attention-deficit/hyperactivity disorder samples and validated the results in zebrafish models. Methods We used 2 relatively large attention-deficit/hyperactivity disorder cohorts, 799 and 776 cases, respectively. WM impairment was assessed using the Rey Complex Figure Test. First, association analyses were conducted at single-variant, gene-based, and gene-set levels. Deeper insights into the biological mechanism were gained from further functional exploration by bioinformatic analyses and zebrafish models. Results Genomic analyses identified and replicated a locus with rs75885813 as the index single nucleotide polymorphism showing significant association with WM defects but not with attention-deficit/hyperactivity disorder. Functional feature exploration found that these single nucleotide polymorphisms may regulate the expression level of RBFOX1 through chromatin interaction. Further pathway enrichment analysis of potential associated single nucleotide polymorphisms revealed the involvement of posttranscription regulation that affects messenger RNA stability and/or alternative splicing. Zebrafish with functionally knocked down or genome-edited rbfox1 exhibited WM impairment but no hyperactivity. Transcriptome profiling of rbfox1-defective zebrafish indicated that alternative exon usages of snap25a might partially lead to reduced WM learning of larval zebrafish. Conclusions The locus with rs75885813 in RBFOX1 was identified as associated with WM. Rbfox1 regulates synaptic and long-term potentiation–related gene snap25a to adjust WM at the posttranscriptional level. Keywords: Genetic basis, RBFOX1, Transcription regulation, Transcriptome profile, Working memory, Zebrafish __________________________________________________________________ Working memory (WM) is the ability to maintain and manipulate information in the brain/central nervous system to guide goal-directed behaviors, requiring gene expression control. WM impairments are a common trait of many psychiatric disorders, especially attention-deficit/hyperactivity disorder (ADHD) ([45]1). ADHD has an estimated prevalence of 5.9% in school-age children; with core symptoms of inattention, hyperactivity, and impulsivity, children with ADHD usually exhibit impairments in WM ([46]2). WM defects are the leading cause of academic failure in patients with ADHD ([47]3), yet the underlying genetic and neurobiological bases are still not fully understood. To date, studies of the etiology of ADHD have identified only a few variants of candidate genes, largely due to the heterogeneous nature of clinical phenotypes ([48]4). Endophenotypes, on the other hand, may play important roles in unveiling the psychopathological processes of most psychiatric disorders. In this sense, WM defect is a reliable and promising endophenotype that could combine various pathways, including neurotransmission systems, ion channels, transcription regulators, and neurodevelopmental genes ([49]5, [50]6, [51]7, [52]8, [53]9). In the neurotransmission system, the genes of presynaptic components, such as calcineurin, DTNBP1, dysbindin, GAD1, and SNAP25 ([54]10, [55]11, [56]12, [57]13, [58]14, [59]15); postsynaptic proteins, such as AKT/AKT1, GRIN1, GRIN2B, NOS1, and NRG1 ([60]16, [61]17, [62]18, [63]19, [64]20); their interacting factors, such as DAO and OLFM3 ([65]21, [66]22, [67]23); and some transcriptional regulators, such as PAWR, TBX1, ZNF804A, DGCR8, and CYFIP1 ([68]24,[69]25), are all associated with WM. Such a research strategy to pick up candidate genes related to WM defects has shed light on the understanding that defects in the neurotransmission system contribute significantly to WM deficiency. Visual WM is crucial in processing visual information. Its deficiencies are linked to general dysfunction in cognition ([70]26) and presented as one of the important symptoms in several psychiatric disorders ([71]27,[72]28). Zebrafish, a kind of vertebrate, has high similarity in physiological structure, growth, and development process with human beings and has shown high conservatism in evolution, making it a hot model organism in the field of biomedical research in recent years ([73]29,[74]30). Although the central nervous system of zebrafish is different from that of mammals, several key brain regions of zebrafish are homologous to that of mammals ([75]31). The behavioral patterns of zebrafish and mammals are also quite similar. Zebrafish can show high-level behavior and neural integration, including memory, conditioned reflex, and social behavior ([76]32,[77]33). In this study, we performed genomic analyses of visuospatial WM in a relatively large number of children with ADHD and validated the genetic variants significantly associated with WM abnormality. We further tested the affected gene in zebrafish, taking advantage of this highly tractable model system ([78]29,[79]30). Methods and Materials Subjects: Discovery and Replication Cohorts Two consecutive ADHD samples were recruited from our child psychiatric clinics, which included 1040 cases in the discovery cohort and 1192 cases in the replication cohort. Both cohorts were medication free and between 6 and 16 years of age. All cases met DSM-IV ADHD diagnostic criteria based on a semistructured interview by senior child and adolescent psychiatrists using the Clinical Diagnostic Interview Scale. The other inclusion criteria were as follows: Full Scale IQ ≥ 70 and both biological parents were of Han descent. Those who had comorbidities with major neurological or psychiatric disorders, such as schizophrenia, bipolar disorder, major depressive disorder, pervasive development disorder, and epilepsy, were excluded. Individuals without WM measures (n = 200 and 342, respectively), age (n = 0 and 2, respectively), or IQ (n = 27 and 37, respectively) were excluded from genetic analyses (see the pipeline in [80]Figure S11). Majority of the cases were medication naïve. If the child had been medicated, a washout period for at least 1 month was necessary before the recruitment. This work was approved by our Institutional Review board. Written informed consent was obtained from parents. Visuospatial WM Task: Rey Complex Figure Test Delayed Component The subject was instructed to observe a complex figure designed by Rey ([81]34) and then was asked to draw the figure from memory onto a blank sheet of paper after 30 seconds. After a 20-minute delay, the subject was asked to recall and draw the figure from memory again. The test was scored according to the structure (0–6) and detail accuracy (0–36). The delayed structure and detail scores were the primary variables of interest in this study for their higher cognitive load. We conducted principal component analyses to extract common features from delayed structure and detail accuracy to produce an index of Rey Complex Figure Test Delayed Component (REYD), which assesses delayed WM. Genomic and Bioinformatics Analyses of WM in Patients With ADHD DNA Extraction and Genotyping Genomic DNA was extracted directly from peripheral blood sample of each subject. For the first cohort, all participants were genotyped in our ADHD genome-wide association study project using Affymetrix 6.0 array (Affymetrix, Inc.). For the second cohort, Illumina Infinium PsychArray (Illumina, Inc.) was used. Both were genotyped at CapitalBio Ltd. Genotypes were called by GENOME STUDIO calling algorithm with the human reference genome (hg19). The same quality control steps were performed. Individuals with per-individual autosomal heterozygosity standard deviation larger than the mean, sex inconsistent with site reports, a per-individual call rate < 95%, and the lower call rate in a pair of individuals with proportion identity by descent PI_HAT > 0.185 were excluded ([82]35). Furthermore, relatedness between the 2 cohorts was checked (PI_HAT < 0.05). Then, the variants were filtered based on per–single nucleotide polymorphism (SNP) call rate < 95%, deviation from Hardy-Weinberg equilibrium with p < .001, or a minor allele frequency < 1% ([83]35). After quality control, 1026 samples with 644,166 SNPs and 1147 samples with 284,176 SNPs remained for the 2 cohorts, respectively (pipeline shown in [84]Figure S11). We imputed nongenotyped SNPs of the 2 cohorts using IMPUTE2, with 2186 phased haplotypes from the full 1000 Genomes Project Integrated Phase 3 Release ([85]36) as the reference panel. We removed imputed SNPs with a squared correlation with the true genotypes r^2 < 0.9 or with minor allele frequency < 0.01. Finally, 6,552,994 and 5,468,003 SNPs were included for the 2 cohorts, respectively, after imputation. Genetic Single-Marker Analyses We performed association analyses on a single SNP for REYD using PLINK version 1.9 ([86]37). By multidimensional scaling, no substantial population stratification was found. As WM is correlated with age, sex, and IQ ([87]38,[88]39), we included them as covariates. Linear regression models were used with age, sex, IQ, and the top 10 eigenvectors from the genetic principal component analysis as covariates for pruned SNPs (r^2 = 0.2) by EIGENSOFT 4.2 ([89]40). When we performed phenotypic association analyses for ADHD (i.e., ADHD diagnosis and inattention, impulsivity, and hyperactivity measures based on Clinical Diagnostic Interview Scale), all covariates mentioned above were included except IQ. Bonferroni correction was used, and p < 5 × 10^−8 was regarded as whole-genome–wide significance. Significant loci were tested in the replication sample (corrected p [p[CORR]] < .05). Meta-analyses were also implemented in PLINK ([90]37). All reported p values are 2-sided. Gene-Based and Pathway Enrichment Analyses We further conducted gene-based analyses for both REYD and ADHD using MAGMA ([91]41). Then, we performed pathway enrichment analyses using SNPs with p < 1 × 10^−4 from the meta-analyses as implemented in MAGMA ([92]41) with a 35-kb upstream and 10-kb downstream window around genes as the default setting. A total of 10,185 Gene Ontology (GO) and 186 Kyoto Encyclopedia of Genes and Genomes gene sets [obtained from MSigDB ([93]42)] were included. Competitive p values were computed and interpreted. MAGMA’s built-in Bonferroni correction was used for multiple testing corrections. Polygenic Risk Score Then, we performed polygenic analyses to detect a shared genetic basis between WM and ADHD symptoms. Based on the association analysis on the REYD in the discovery sample, we used PRSice ([94]43) to perform the polygenic risk score analyses to select the most precise threshold for the p value that predicted ADHD symptoms in both discovery and validation cohorts, respectively, with step increased p-value thresholds (.000001, .00001, .0001, .001, .01, .02, .03, .04, .05, .1, .2, .3, .4, .5, and 1). An empirical multiple testing correction implemented in PRSice was applied, which is based on a permutation procedure. The significance of the regression results was corrected by a permutation test with 10,000 replicates, and α = 0.001 as suggested was used ([95]43). We also explored whether genetic components of glutamate receptor or long-term potentiation (LTP)–related pathways were shared between WM and ADHD symptoms. We included 236 and 75 genes that belong to glutamate receptor and LTP-related pathways, respectively, according to MSigDB ([96]42) and performed the same polygenic risk score analyses as mentioned above. Regulatory Feature Analyses and Network Construction We obtained the regulatory features of the significant SNPs from rSNPBase version 3.0 ([97]44) and HaploReg version 4.1 ([98]45), two user-friendly graphic interface web tools that integrate comprehensive information related to genomic regulation. We also searched the expression quantitative loci in the Genotype-Tissue Expression portal ([99]46). The expression plot was generated using the Human Protein Atlas. Using STRING ([100]47), the interacting genes network was constructed for RBFOX1 and those enriched genes from pathway analyses. Functional Analyses of rbfox1 in Zebrafish Zebrafish Husbandry and Care Zebrafish were raised at 28 °C with a density of 8 to 10 fish/L and experimented with the established standards ([101]48). The wild-type zebrafish used in this study was Oregon AB strain/line. All zebrafish experiments were conducted according to the guidelines approved by the Institutional Animal Care and Use Committee of Southern University of Science and Technology. Zebrafish Loss-of-Function Studies Loss-of-function (knockdown or knockout) experiments conducted in this study were achieved by morpholino antisense oligo technique (Gene Tools Inc.), 4 guide RNAs (gRNAs)/Cas9 technique, or Cas9 ribonucleoproteins technique. To obtain mutants of zebrafish rbfox1 to evaluate if the gene affects visuospatial memory, we used synthetic crRNA:tracrRNA (crispr RNA:transactivating crispr RNA) duplex gRNAs and Cas9 ([102]49). Two target sequences in RNA recognition motif of Rbfox1 were chosen, and each created indels (found in each copy of the target locus of rbfox1) in more than 80% of injected embryos ([103]Figure S8). Because most injected F0 embryos could be treated as true null mutants ([104]49), we used 5 days post fertilization (dpf)–injected larval zebrafish and morpholino antisense oligo-injected morphant to conduct behavior tests. In the case of snap25a, 4 gRNAs ([105]Table S4) ([106]http://crispor.tefor.net/) were coinjected with Cas9 nuclease, and the resulting F0 fish larvae were called knockout fish or snap25a 4-gRNA fish. Genotyping Genomic DNA was extracted from 24 hours post fertilization zebrafish embryos using TIANamp Genomic DNA Kit (Tiangen). The aimed region(s) of the gene locus was amplified by gene-specific polymerase chain reaction (PCR) primers, and the expected DNA fragment was purified by QIAquick Gel Extraction Kit (Qiagen). The purified PCR product(s) was cloned for sequencing purpose (Sangon Biotech). The cloning primers sequence are listed in [107]Table S5. Quantitative Real-Time PCR Total RNA was extracted from 50 wild-type or mutant larvae at 5 dpf after fixing them in TRIzol reagent (Invitrogen). A total of 1 μg RNA was reverse-transcribed to the first strand of complementary DNA with the random primer using a complementary DNA synthesis kit (Promega). Quantitative real-time PCR primers for snap25a and β-actin gene are listed [108]Table S6. The quantitative reverse transcription PCR was performed in an ABI 7500 Real-Time PCR instrument (Applied Biosystems) with the SYBR green detection system, and results were normalized with β-actin expression using ΔΔCt method. Transcriptome Sequencing Analysis Total RNA was separately prepared from 50 larvae (5 dpf) of wild-type or rbfox1 morphant and sequenced through the DNBSEQ platform, with the National Center for Biotechnology Information accession number GCF_000002035.6_GRCz11 as the zebrafish genome reference. The sequencing data were filtered with SOAPnuke (version 1.5.2), and then the clean reads were mapped to genome using HISAT2 (version 2.0.4). After aligning the clean reads to genes via Bowtie2 (version 2.2.5), the Beijing Genomics Institute created a database to include all annotated coding transcripts with actual number of reads/counts calculated by RSEM (version 1.2.12). Based on the hypergeometric test and corrected by Q value with a rigorous threshold (Q ≤ 0.05), respectively, the GO and Kyoto Encyclopedia of Genes and Genomes enrichment analysis of differential gene expressions and significant levels of terms and pathways were obtained. Zebrafish Behavioral Tests: WM WM can be represented by habituation in zebrafish. Short-term habituation could be analyzed through repeated acoustic or visual stimulations. Free swimming 5-dpf zebrafish larvae were dispersed in a 96-well plate to sit in the observation chamber (DanioVision; Noldus). Response velocity to acoustic/vibrational or dark flash stimuli was detected and calculated by using EthoVision XT13 video-tracking software (Noldus). For the acoustic habituation assay, 10 stimuli with an acoustic intensity of 90 dB were delivered with 1-second interstimulus interval (ISI) as baseline, followed by 20 stimuli with 1-second ISI. Acoustic habituation was indicated as % acoustic habituation = [1 − (velocity of stimuli 21–30)/(response velocity of baseline)] × 100. For the visual habituation assay, a 4-block training protocol was performed without any break between blocks. Each block consisted of 120 dark flashes with a 15-second ISI, and each dark flash lasted for 1 second. Habituation was indicated as % visual habituation = (1 – block 4/block 1) × 100 ([109]Figure S7). Statistical Analyses of Zebrafish Experiments GraphPad Prism version 7.00 was used for statistical analysis of zebrafish experiments. For 2-group comparisons, one-tailed t test with 95% CI was used. For comparisons with 3 or more groups, one-way analysis of variance with 95% CI was used. Results Genome-wide Association of WM and the Relationship With ADHD Symptoms at Single-Variant and Polygenic Levels We first performed a genome-wide analysis of visuospatial WM in children with ADHD. A total of 799 and 776 cases were ultimately included in the discovery and replication stages, respectively (multidimensional scaling plots are available in the [110]Supplement). The demographic, IQ, and cognitive phenotype data for the discovery and replication samples are presented in [111]Table S1. We used delayed component (REYD), which was the principal component axis score accounting for 94.1% and 91.1% variance in the discovery and validation cohorts, respectively, in the following analyses. In the discovery stage, we identified 5 significant SNPs in whole-genome–wide association analysis, among which rs75885813 was ranked on top (p[CORR] = 3.83 × 10^−9). The other 4 SNPs were in high linkage disequilibrium with rs75885813 (r^2 > 0.9) ([112]Table 1). Quantile-quantile plot for SNP associations is presented in [113]Figure S2. All 5 SNPs were located on 16p13.3 (Chr16:7120001-7160000, hg19) within the RBFOX1 gene ([114]Figure 1). They reached statistical significance in both the replication stage (p[CORR] < .05) and meta-analyses (p[CORR] < 5 × 10^−9). However, associations between the RBFOX1 SNPs and core symptoms of ADHD (i.e., inattention, hyperactivity, impulsivity, and overall symptoms) were not statistically significant (p[CORR] > .05, data available upon request), nor was the association for PRS weighted by WM for patients with ADHD (p[CORR] > .05) ([115]Table S2). Table 1. Significant SNPs Identified From Genome-wide Analysis of REYD in Discovery and Replication Samples and Meta-analyses CHR SNP Pos A1 LD With rs75885813 β __________________________________________________________________ SE __________________________________________________________________ Coefficient t-Statistic __________________________________________________________________ p __________________________________________________________________ Discovery Replication Discovery Replication Discovery Replication Discovery Replication Meta-analysis 16 rs75885813 7,141,263 T – −0.669 −0.276 0.112 0.130 −5.96 −2.12 3.83 × 10^−9 0.034 3.67 × 10^−9 16 rs114891671 7,142,203 C 0.98 −0.669 −0.352 0.113 0.135 −5.90 −2.60 5.56 × 10^−9 0.0094 6.04 × 10^−10 16 rs115712598 7,142,236 T 0.98 −0.669 −0.352 0.113 0.135 −5.90 −2.60 5.56 × 10^−9 0.0094 6.04 × 10^−10 16 rs116183707 7,142,287 A 0.98 −0.669 −0.352 0.113 0.135 −5.90 −2.60 5.56 × 10^−9 0.0094 6.04 × 10^−10 16 rs7205995 7,140,331 G 0.98 −0.667 −0.35 0.113 0.135 −5.88 −2.59 6.05 × 10^−9 0.0097 6.69 × 10^−10 [116]Open in a new tab Five whole-genome–significant SNPs that are in high LD were identified (r^2 > 0.8). CHR, chromosome; LD, linkage disequilibrium; Pos, position; REYD, Rey Complex Figure Test Delayed Component; SNP, single nucleotide polymorphism. Figure 1. [117]Figure 1 [118]Open in a new tab Regional plot of significant locus for genome-wide association of Rey Complex Figure Test Delayed Component. The gray horizontal line represents the threshold for genome-wide significant association (p = 5 × 10^−8). rs75885813, corrected p = 3.83 × 10^−9, was the top-ranked SNP. The other 4 significant SNPs were in high linkage disequilibrium with rs75885813 (r^2 > 0.8). All significant SNPs were located within RBFOX1. SNP, single nucleotide polymorphism. Associated Variants May Affect Transcriptional Regulation Given the strong association between rs75885813 (i.e., an intron variant of RBFOX1) and REYD, we further examined the regulatory features of this variant, as well as its interaction networks and coexpression genes. The chromatin state analyses showed that rs75885813 is in functional regulatory regions in various tissues indicated by the enhancer-specific H3K4me1 and promoter-specific H3K4m3 markers. It also alters the regulatory motifs of some transcription factors in the brain ([119]Table S3). The Hi-C data suggested that a chromatin loop can form between the genomic region (16:7120001-7160000, hg19) where rs75885813 is located and the promoter of RBFOX1 (16:6040001-6080000, hg19), also suggesting a possible regulatory effect ([120]Figure S3) on RBFOX1 transcriptional regulation. We further constructed a network to include its interacting partners and coexpressing genes. As shown in [121]Figure S4, the functional protein interaction analyses indicated that most of them (8 of 11) are RNA-binding proteins and may play roles in RNA alternative splicing events. Gene-Based and Pathway Enrichment Analyses No gene achieved significance after multiple corrections for either REYD or ADHD. We investigated subthreshold variants from the REYD association analysis using a pathway enrichment test. After removing gene sets that contained <2 genes defined in the GO analyses, 773 gene sets were included. Sixteen GO pathways reached significance (p[CORR] < .05) ([122]Figure 2 and [123]Table 2). All of these pathways are involved in posttranscription regulation, among which the mRNA metabolic process was ranked on top. The relationships of these significantly enriched pathways are presented in [124]Figure S5. The extended interaction gene network showed the potentially interacted genes, many of which are overlapped with risk genes for psychiatric disorders ([125]Figure S6). Figure 2. [126]Figure 2 [127]Open in a new tab Bubble diagram of the pathway analysis of potentially associated variants for Rey Complex Figure Test Delayed Component. The p value increased from top to bottom. All the pathways displayed in the figure were significant (corrected p < .05). The color represents p value, while the size represents the number of associated genes enriched in the particular pathway. mRNA, messenger RNA. Table 2. Enriched Biological Process by Potential Associated Variants With Working Memory Full Name No. of Genes/Total[128]^a Gene List β Standardized β SE p GO:0016071 GOBP mRNA Metabolic Process 5/879 RBFOX1, SMG6, TNFSF13, FXR2, RAVER1 0.697 0.144 0.0995 6.963 × 10^−10 GO:1903311 GOBP Regulation of mRNA Metabolic Process 3/334 RBFOX1, TNFSF13, FXR2 1.511 0.244 0.235 7.520 × 10^−9 GO:0098791 GOCC Golgi Apparatus Subcompartment 2/887 RBFOX1, CNGB1 1.454 0.193 0.230 1.189 × 10^−8 GO:0003729 mRNA Binding 3/538 RBFOX1, EIF4A1, FXR2 1.416 0.229 0.231 2.538 × 10^−6 GO:0006397 GOBP mRNA Processing 3/543 RBFOX1, FXR2, RAVER1 1.350 0.218 0.222 2.917 × 10^−8 GO:0000375 GOBP RNA Splicing via Transesterification Reactions 3/384 RBFOX1, FXR2, RAVER1 1.350 0.218 0.222 2.917 × 10^−8 GO:0051236 GOBP Establishment of RNA Localization 3/200 NUP205, RBFOX1, SMG6 0.590 0.0952 0.0979 3.882 × 10^−8 GO:0006403 GOBP RNA Localization 3/233 NUP205, RBFOX1, SMG6 0.590 0.0952 0.0979 3.882 × 10^−8 GO:0000381 GOBP Regulation of Alternative mRNA Splicing via Spliceosome 2/57 RBFOX1, FXR2 1.370 0.181 0.228 4.009 × 10^−8 GO:0050684 GOBP Regulation of mRNA Processing 2/139 RBFOX1, FXR2 1.370 0.181 0.228 4.009 × 10^−8 GO:0048024 GOBP Regulation of mRNA Splicing via Spliceosome 2/100 RBFOX1, FXR2 1.370 0.181 0.228 4.009 × 10^−8 GO:0043484 GOBP Regulation of RNA Splicing 2/144 RBFOX1, FXR2 1.370 0.181 0.228 4.009 × 10^−8 GO:0035770 GOCC Ribonucleoprotein Granule 2/244 RBFOX1, FXR2 1.370 0.181 0.228 4.009 × 10^−8 GO:0015931 GOBP Nucleobase Containing Compound Transport 4/252 NUP205, SLC28A1, RBFOX1, SMG6 0.511 0.0948 0.0946 4.537 × 10^−7 GO:0003723 GOMF RNA Binding 9/1938 RPN1, C7orf50, RALYL, RPP25L, RBFOX1, SMG6, EIF4A1, FXR2, RAVER1 0.491 0.134 0.0940 8.876 × 10^−7 GO:0008380 GOBP RNA Splicing 5/478 C2orf49, RBFOX1, USB1, FXR2, RAVER1 1.091 0.225 0.219 2.360 × 10^−6 [129]Open in a new tab The pathway enrichment analyses were conducted with single nucleotide polymorphisms with p < 1 × 10^−4 from the association analyses as implemented in MAGMA. GOBP, Gene Ontology biological processes; GOCC, Gene Ontology cellular components; GOMF, Gene Ontology molecular functions; mRNA, messenger RNA. ^a The number of genes enriched in the pathway and the total number of genes of the pathway. Functional Analyses of rbfox1 in Zebrafish To test our hypothesis that a downregulation of Rbfox1 might lead to WM defect, we first established a measurement to monitor habituation learning behaviors in zebrafish larvae. We found that during acoustic/tapping stimulation with a total of 30 taps, larval zebrafish gradually reduced the extent of their startle responses ([130]Figure S7A, B). After exposure to a massed/continued dark flash-training period, visual habituation was also evident ([131]Figure S7C, D). Mutant-like larval zebrafish (5 dpf), coinjected with chemically modified crRNAs/tracrRNAs and Cas9 protein against 2 sites of rbfox1 at the 1-cell stage in wild-type fertilized eggs ([132]Figure S8), exhibited a decrease in habituation ([133]Figure 3A–D). The knockdown of rbfox1 function with morpholino antisense oligos that blocked rbfox1 pre-mRNA splicing resulted in similar habituation phenotypes ([134]Figure 3E, F), suggesting a relationship between Rbfox1 and the WM defect in zebrafish. Figure 3. [135]Figure 3 [136]Open in a new tab rbfox1 is required for visual and acoustic habituation in larval zebrafish. (A) Schematic diagram of visual habituation. (B) Mean percentage of habituation during the visual habituation phase in the rbfox1 mutant. (C) Schematic diagram of acoustic habituation. (D) Mean percentage of habituation during the acoustic habituation phase in the rbfox1 mutant. (E) Mean percentage of habituation during the visual habituation phase in the rbfox1 MO-injected fish. (F) Mean percentage of habituation during the acoustic habituation phase in the rbfox1 MO-injected fish. The data are expressed as mean ± standard error. n = 48 larval fish per group. ∗p < .05 (paired t test), ∗∗p < .01. MO, morpholino antisense oligo; mut, mutant; WT, wild-type. Posttranscriptional Regulation of WM Genes To select zebrafish WM-related genes in an unbiased way, we profiled the transcriptome of rbfox1 splicing morpholino antisense oligo-injected larvae at 5 dpf. The GO and Kyoto Encyclopedia of Genes and Genomes enrichment analyses indicated that the most significantly changed pathways included the LTP, which was tightly related to memory formation, and transient receptor potential channels, which had been reported to define hippocampal synaptic transmission and WM ([137]Figure 4A). The most affected cellular components included pre- and postsynapse connections ([138]Figure 4B). Gene expression of proteins in postsynaptic density, the presynaptic membrane, and the ionotropic glutamate receptor complex were all significantly changed. For instance, we found that snap25a and snap25b were the most abundantly expressed genes encoding a synaptosome-associated protein in 5-dpf larvae ([139]Figure S9). Figure 4. [140]Figure 4 [141]Open in a new tab Zebrafish rbfox1 regulates both presynaptic and postsynaptic functions. (A) Kyoto Encyclopedia of Genes and Genomes pathway enrichment of differentially expressed genes in rbfox1 morphant. (B) Gene Ontology enrichment of differentially expressed genes in rbfox1 morphant. All zebrafish larvae were collected at 5 days post fertilization, and then whole transcriptome deep sequencing was performed. HIF, hypoxia-inducible factor; GnRH, gonadotropin-releasing hormone; Num, number; TRP, transient receptor potential. Among well-studied neurotransmission-related genes, functional protein association analyses revealed that SNAP-25 is within a few closely related RNA processing components networked with RBFOX1 ([142]Figure S4). Zebrafish Snap25a might be the most abundant synaptosome-associated protein in 5-dpf larval fish responsible for presynaptic neurotransmitter release ([143]Figure S9A). Zebrafish snap25a is a homolog of mouse Snap25, which is known to have 2 alternative splicing variants: Snap25a and Snap25b that correspond to 2 zebrafish splicing variants, snap25a-202 and snap25a-201 transcripts ([144]Figure S9D), respectively, during 24 to 96 hours post fertilization. Only 11 amino acids are different between Snap25a-201 and Snap25a-202, and these amino acids are all due to the alternative usage of exon-5. Apparently, snap25a-201 is predominantly expressed during embryogenesis and early larval stage ([145]Figure S9C). In rbfox1 mutant, snap25a-201 transcripts were significantly decreased (p < .01) while snap25a-202 expression was significantly increased (p < .001) ([146]Figure 5B). An rbfox1 binding motif, (U)GCAUG, is found in exon 5 of snap25a-202, and the overall reduced snap25a expression level ([147]Figure S9B) also indicates that an imbalanced expression of two splicing variants most likely changes zebrafish visuospatial memory or acoustic/vibration-mediated memory. Figure 5. [148]Figure 5 [149]Open in a new tab Posttranscriptional regulation of snap25a by Rbfox1. (A) Acoustic and visual habituation defects in snap25a-4sg mutant, tests performed on 5 days post fertilization, n ≥ 48. (B) Changes of snap25a isoform ratios (exon5 alternative inclusion) in rbfox1 mutant. The data are expressed as mean ± standard error. n = 48 larval fish per group. ∗∗∗p < .001, ∗∗p < .01, ∗p < .05 (unpaired t test). 4gRNA, 4 guide RNA; mut, mutant; WT, wild-type. To obtain zebrafish snap25a mutants, we used another gene knockout method that produces null F0 zebrafish with high probability ([150]50). Using a mix of 4 single-guide RNAs against snap25a and Cas9 to inject the yolk of fertilized eggs, we evaluated learned memory of snap25a mutant at 5 dpf and discovered significantly reduced acoustic/vibration and visual habituation ([151]Figure 5A), similar to those found in rbfox1 mutant and morphant that were responsible for the impaired short-term memory. Discussion The past few decades have witnessed studies of many candidate genes related to the genetic vulnerability of WM impairments. However, most picked candidate genes were limited to the hypotheses. Our genome-wide single-variant and polygenic analyses of WM defects in children with ADHD identified 5 RBFOX1 intronic variants/SNPs that may affect RBFOX1 levels in the brain. Zebrafish loss-of-function of rbfox1 experiments showed both visual and auditory habituation defects in larval fish. When rbfox1 was knocked down or out, snap25a was found downregulated and mis-spliced. To our knowledge, this is the first study to uncover that RBFOX1/rbfox1 regulates WM posttranscriptionally. Our pathway enrichment analysis has also identified transcription regulation that is linked to WM in accordance with previous candidate gene, gene-set, and polygenic studies ([152]24,[153]25). Some of the regulators are important for the synthesis of proteins required during neurobiological process of WM and mRNA metabolic process. If the expression of more RBFOX1/rbfox1-regulated LTP/synaptic genes is found to be controlled at the pre-RNA splicing and mRNA stability stages, the fine and quick regulation of neural activities such as habituation, startle response, and dependent transcription could be economically and efficiently achieved. RBFOX1/rbfox1 is a pleiotropic gene that has been associated with 7 specific psychiatric disorders ([154]51). The most significant SNP, rs75885813, was previously associated with Alzheimer’s disease-related phenotypes in patients with Alzheimer’s disease whose characteristic traits include WM deficits ([155]52,[156]53). The association between RBFOX1 rs7193263 and major depression is also evident in a genome-wide association study (p = 9.73 × 10^−9) ([157]54). RBFOX1 and WM were never directly associated, yet in a genome-wide gene expression study, Rbfox1 was linked to RNA processing after memory retrieval ([158]55). However, in this study, no association was found between hyperactivity, an ADHD core symptom, and RBFOX1/rbfox1 or snap25a, either in the population genetic analyses or in animal experiments ([159]Table S2 and [160]Figure S10). We noticed that the association of RBFOX1 was mainly driven by a major depression sample (genome-wide association study p = 9.73 × 10^−9) among 7 psychiatric disorders ([161]54). The p value for the association of ADHD was only .0065 ([162]56). Thus, the effect of RBFOX1 on WM is potentially not confined only to ADHD. This study further implies that RBFOX1 may be the common susceptibility gene of psychiatric disorders via regulating WM. More association studies of RBFOX1-regulated neurotransmission system genes and WM may be needed to strengthen our view. In the functional protein association networks of RBFOX1, SNAP25 attracted our attention. A previous study indicated that SNAP25 is related to WM deficits in patients with ADHD ([163]14). In our study, we found a reduction in snap25a and snap25b expression in zebrafish rbfox1 morphant, which showed WM defect but no hyperactivity. The decrease of snap25a expression was more pronounced than that of snap25b. Zebrafish snap25a is homologous to mouse Snap25, which is known to have 2 alternative splicing variants: Snap25a (exon 5a, expressed late, and corresponding to zebrafish snap25a-202) and Snap25b (exon 5b, expressed early, and corresponding to zebrafish snap25a-201). Mouse Snap25a and Snap25b are expressed during embryonic and early postnatal development, respectively. Snap25 differentially affects interactions with other SNAREs (SNAP receptors) and SNARE-interacting proteins ([164]57,[165]58). The Rbfox1^−/− mouse brain exhibited normal mRNA level, but decreased Snap25b transcripts and increased Snap25a transcripts ([166]59), consistent with our findings of a decreased snap25a-201/25a-202 ratio ([167]Figure 5B). It is worth exploring whether reversing the snap25a-201/25a-202 ratio or simply increasing snap25a-201 mRNA can rescue the WM defect in rbfox1 mutants in future studies. However, there were some limitations in this study. First, due to the limited sample sizes, this study is underpowered to detect other WM-related variants and unveil its genetic structure more comprehensively. Enlarging the sample size and cooperation with other groups could strengthen the statistical power. In summary, these findings have revealed possible causal pathways of postranscription regulators that trigger WM deficits and alternative splicing events, mediated by Rbfox1, control SNARE, and LTP genes to affect WM. Acknowledgments and Disclosures