Plausibility of the zebrafish embryos/larvae as an alternative animal model for autism: A comparison study of transcriptome changes

Autism spectrum disorder (ASD) is a serious neurodevelopmental disorder characterized by impaired or abnormal social interaction and communication and by restricted and repetitive behaviour. ASD is highly prevalent in Asia, Europe, and the United States, and the frequency of ASD is growing each year. Recent epidemiological studies have indicated that ASD may be caused or triggered by exposure to chemicals in the environment, such as those in the air or water. Thus, toxicological studies are needed to examine chemicals that might be implicated. However, the experimental efficiency of existing experimental models is limited, and many models represent challenges in terms of animal welfare. Thus, alternative ASD animal models are necessary. To address this, we examined the efficacy of the zebrafish embryo/larva as an alternative model of ASD. Specifically, we exposed zebrafish to valproic acid (0, 12.5, 25, 50, or 100 μM), which is a chemical known to induce autism-like effects. We then analysed subsequent developmental, behavioural, and transcriptomic changes. We found that 100 μM and 50 μM valproic acid decreased the hatching rate and locomotor activity of zebrafish embryos/larvae. Transcriptomic analysis revealed significant alterations in a number of genes associated with autism, such as adsl, mbd5, shank3, and tsc1b. Additionally, we found changes in gene ontology that were also reported in previous studies. Our findings indicate that zebrafish embryos/larvae and humans with ASD might have common physiological pathways, indicating that this animal model may represent an alternative tool for examining the causes of and potential treatments for this illness.


Introduction
Autism spectrum disorder (ASD) is an early onset neuropsychiatric disorder that is characterized by impaired social communication, restricted interests, and stereotyped and repetitive behaviours [1,2]. Recent years have seen a dramatic increase in the prevalence of ASD in Asia, Europe, and the USA. The prevalence of ASD in Sweden (Stockholm cohort study) increased by almost 3.5-fold between 2001 and 2011 [3], and the prevalence in the United States increased from 1 in

Test chemical and exposure experiments
Valproic acid sodium salt (2-propylpentanoic acid sodium, CAS RN. 1069-66-5, purity > 98%) was purchased from Sigma-Aldrich (St. Louis, MO, USA). The test solutions were prepared in E3 medium, 0.292 g NaCl, 0.013 g KCl, 0.044 g CaCl, and 0.081 g MgSO 4 in Millipore water (1 L). The exposure concentrations (0, 12.5, 25, 50, 100 μM of VPA) were determined according to previous studies [31,34]. Normal zebrafish embryos were randomly distributed into test solutions of VPA for 120 h. All exposure experiments were conducted in accordance with protocols approved by the Institutional Animal Care and Use Committee (IACUC) of the Korea Institute of Toxicology (Protocol No. RS16005).

Survival, hatching status, and malformation
During the exposure period, embryo/larva survival, hatching rate, time to hatch, and morphological malformation rate were checked every 24 h until 120 h. Percentages corresponding to mortality, hatching rate, and malformation rate were determined in quadruplicate (N = 4), with six embryos in each sample.

Locomotor activity (movement tracking)
After 120 h exposure to VPA, we monitored the movement of each fish via an automated tracking device, Zebrabox™ (Viewpoint, Lyon, France). We tracked the larvae per treatment with octuplicates (N = 8) including three larvae in each. Tracking occurred during 4 repeated phases of a 3-min light/3-min dark period for a total of 24 min. The video data were then analysed using ZebraLab™ software to calculate the total distance travelled and active duration rate ((total movement duration-inactive duration) / total duration) for each individual larva.

RNA-seq analysis
After 72 h and 120 h exposure, we pooled twenty zebrafish larvae per sample for RNA extraction. Triplicate (N = 3) pooled RNA samples from each treatment group were prepared and kept in RNAlater RNA stabilization solution (Qiagen, Hilden, Germany) for further RNA extraction. Total RNA was isolated using the RNeasy mini kit (Qiagen, Hilden, Germany). We assessed rRNA band integrity using an Agilent RNA 6000 Nano kit (Agilent Technologies, CA, USA). We used samples with an RNA Integrity Number (RIN) greater than 7 for RNA library construction. Briefly, prior to cDNA library construction, we used 1 μg of total RNA and magnetic beads with oligo (dT) to enrich the poly (A) mRNA. Then, the purified mRNA was disrupted into short fragments, and the double-stranded cDNA was immediately synthesized. The cDNA was subjected to end-repair, poly (A) addition, and connected with sequencing adapters using the TruSeq RNA Sample Prep Kit (Illumina, Ca, USA). The suitable fragments, purified via a BluePippin 2% agarose gel cassette (Sage Science, MA, USA), were selected as templates for PCR amplification. The final libraries were quantified using a KAPA library quantification kit (KAPA Biosystems, South Africa), and the quality of the library was evaluated using an Agilent 2100 bioanalyser (Agilent Technologies, CA, USA). The fragments were found to contain between 350 and 450 base pairs. Subsequently, the library was sequenced using an Illumina HiSeq2500 sequencing platform (Illumina, CA, USA). Low-quality reads were filtered according to the following criteria: reads contained more than 10% skipped bases, reads contained more than 40% of bases whose quality scores are less than 20, and reads with average quality scores of each read less than 20. The filtered reads were mapped to the zebrafish reference genome (Ensembl version 86) using the aligner STAR v.2.3.0e [35].

DEG analysis and GO analysis
We measured the gene expression level using Cufflinks v2.1.1 [36] with the gene annotation database from Ensembl version 86. The non-coding gene region was removed with the mask option. To improve the accuracy of measurement, we applied multi-read correction and frag bias correction. The abundance of gene transcripts was measured via FPKM (fragments per kilobase of transcript per million fragments mapped). The FPKM cut-off level was set at 0.1. For differential expression analysis, gene level count data were generated using the HTSeq-count v0.5.4p3 tool [37]. Using calculated read count data, differentially expressed genes (DEGs) were identified using the TCC R package [38]. Genes with a p-value < 0.05 and | log 2 FC(fold change) | > 1 were considered to be differentially expressed. Significantly altered DEGs were compared with the ASD-related gene lists from previous studies [2,6,9,11,39]. To predict the function of the selected genes, the DEGs were subjected to gene ontology (GO) according to three categories: biological process (BP), cellular component (CC), and molecular function (MF). The ontology and annotation files for GO enrichment analysis were downloaded from the gene ontology website (http://www.geneontology.org/). p-values < 0.001 were considered statistically significant [40]. Significantly altered GOs were compared with the ASD-related GO lists from previous studies [6,41]. Concentration-dependent changes in GO enrichment were determined based on the gene lists showing transcriptional changes over 2 (| log 2 FC | > 1) between at least two serial concentrations with consistent positive or negative slopes.

Quantitative PCR analysis
To validate the transcriptions of DEG analysis, quantitative PCR analysis was carried out. Target genes, i.e., adsl (adenylosuccinate lyase), mbd5 (methyl-CpG binding domain protein 5), shank3a (SH3 and multiple ankyrin repeat domains 3a), and tsc1b (tuberous sclerosis complex 1), were chosen based on the commonality with previous studies (Table 1). After 120 h of VPA exposure, total RNA was extracted using a Maxwell 16 total RNA purification kit (Promega Co., Madison, WI, USA), and its concentration and quality were checked with an ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, ED, USA). During the RNA extraction, genomic DNA was selectively removed with the clearing agent that was included in the purification kit. cDNAs were synthesized by reverse transcription using an iScript cDNA synthesis kit (BioRad, Hercules, CA, USA). The cDNA concentrations were also measured using an ND-1000 spectrophotometer. The samples were then diluted with purified water to 300 ng/μL. Target genes were amplified using an ABI 7300 RT-PCR system (Applied Biosystems, Foster City, CA, USA). The PCR reaction mixture (total 20 μ/L) contained 10 μL of ABI SYBR Green Master mix (ABI, USA), 1.8 μL of each primer (10 pmol), 4.4 μL of purified PCR-grade water, and 2 μL of cDNA sample. The thermal cycle profile was as follows: pre-incubation at 95˚C for 10 min, 40 cycles of amplification at 95˚C for 10 s, 60˚C for 20 s, and 72˚C for 20 s. The comparative Ct methods were employed to calculate the relative level of gene transcription [42]. The sequence of primers for the target genes and reference gene (β-actin) for zebrafish are shown in S1 Table.

Statistical analysis
The normality and homogeneity of each variance were analysed using the Shapiro-Wilk test and Levene's test. Then, we conducted one-way analyses of variance (ANOVA) with Dunnett's test or T3 tests using SPSS 12 for Windows1 (SPSS, Chicago, IL, USA) to determine significant differences between the control and VPA exposure groups. For non-parametric data, we used the Kruskal-Wallis test. Differences with p-values less than 0.05 were considered significant. For trend analysis, linear regression analysis was carried out, and the statistical significance of a linear trend in the slope was determined (p for trend).

Effects of VPA on survival, hatching status, and development
The survival of zebrafish larvae was unchanged for exposure up to 100 μM of VPA (Fig 1(A)). However, the hatching rate significantly decreased at 100 μM of VPA (Fig 1(B)) to only 29%. Time to hatch was significantly delayed at 50 μM of VPA (Fig 1(C)). We observed no significant developmental changes for exposure up to 50 μM of VPA (Fig 1(D)).

Effects of VPA on larval behaviour
VPA exposure at 50 μM caused a significant decrease in the total distance moved compared with control values (Fig 2(A)). At 50 μM VPA, the active duration ratio during the dark phase was significantly reduced to 20% of that in the control group (Fig 2(B)). However, at 25 μM VPA, the total distance moved and active duration ratio during the dark phase increased by 30% and 25%, respectively, although this was not statistically significant (Fig 2(A) and 2(B)).

Illumina sequencing and mapping of RNA sequencing reads onto the zebrafish genome
High-throughput sequencing generated approximately 76.5-90.8 million (M) pairs of raw reads for each group (S2 Table). From each group, 75.5-89.8 M pairs of clean reads were extracted. Then, these high-quality reads were mapped to the reference zebrafish genome over 82.3%-88.3%. The uniquely mapped reads ranged from 78.8% to 85.3%.

Analysis of differentially expressed zebrafish genes following VPA exposure
The number of DEGs ranged from 921 (12.5 μM at 72 h) to 1468 (25 μM at 120 h). The number of DEGs at 25 μM and 50 μM VPA were greater than those at 12.5 μM VPA at 72 h exposure (S1 Fig). At 120 h, the number of upregulated genes was greater than that of downregulated genes. However, there were more downregulated genes than upregulated genes at 72 h for all VPA conditions except 25 μM. Between 26.6% and 36.4% of the DEGs (> 2-fold change) in the present study matched ASD-related genes reported in previous studies [2,6,9,11] (Fig 3(A)). Among them, 24 DEGs were significant (p < 0.05), as shown in Table 1. The transcription of genes indicated in previous studies to be ASD-related (e.g., adsl, eif4a1 (eukaryotic translation initiation factor 4a1), mbd5, nrxn2 (neurexin2), shank3, and tsc1, etc.) was significantly altered. The detailed data for DEGs that were significantly differentially expressed or expressed over 2-fold (| log2FC | > 1) in previous studies are shown in S3-S7 Tables.
The transcription of genes that had been commonly suggested as ASD-related genes by more than two previous studies and our study were significantly changed (Fig 4). Significant upregulation (adsl and mbd5) and downregulation (shank3a and tsb1b) were observed after 120 h VPA exposure in the qPCR analysis. Significant trends (p for trend) were also shown in the transcription of adsl, mbd5a, and shank3a.

Analysis of GO enrichment in zebrafish after VPA exposure
In this study, we identified 44 GO enrichments that had been classified as ASD-related in previous studies [6,41] (Table 2, S8 and S9 Tables). Half of the GOs were included in the CC category, and 13 and 9 of them were included in the MF and BP categories, respectively ( Table 2). Significantly affected GO enrichments (p < 0.001) in this study were reported in 5.8% to 29.4% of previous human and mouse studies (Fig 3(B)). After VPA exposure, we found significant changes in GO enrichments related to transport (GO: 0006810), intracellular function (GO: 0005622), cytoplasm (GO: 0005737), mitochondria (GO: 0005739), cytosol (GO: 0005829), cytoskeleton (GO: 0005856), and cell projects (GO: 0042995), which are commonly found in the mouse cortex and hippocampus [41].
We found that VPA concentration dependently altered the number of ASD-related GOs (Table 3). A total of 59 GOs were identified as ASD-related, including 25 GOs for the BP, 20 GOs for the CC, and 14 GOs for the MF. Synaptic transmission (GO: 0007268), transmission of nerve impulses (GO: 0019226), multicellular organismal signalling (GO: 0035637), and neurological system processes (GO: 0050877) were included in the BP category. Synapse-related GOs, i.e., synapses (GO: 0045202) and synapse parts (GO: 0044456), were in the CC category, and channel activity-related GOs, such as ion channel activity (GO: 0005216), cation channel activity (GO: 0005261), or ligand-gated ion channel activity (GO: 0015276), were found in the MF category.

Plausibility of zebrafish embryos/larvae as an ASD model based on transcriptome changes
We found a common transcriptional profile (DEGs and GOs) that was also found in previous ASD studies using human populations and mouse models. Thus, our transcriptome profile comparison data support the potential of zebrafish as an alternative animal model for ASD research. The zebrafish model has been proposed for ASD research based on behaviour phenotype [25,33], and several ASD-related genes have been detected in the zebrafish [30]. However, this is the first comparison study to use transcriptome profiles. Our findings may be useful in the development of an alternative and rapid experimental model for screening ASD-related effects induced by environmental chemical exposure.
We observed 24 DEGs that had previously been suggested to be ASD-related genes in zebrafish embryos/larvae after VPA exposure (Table 1). A total of 136 ASD-related genes, presented in S3-S7 Tables, exhibited significant (p < 0.05) or a 2-fold increase in expression (| log2FC | > 1). These results indicate that zebrafish might, at least partly, have common genes and physiological pathways for ASD. The DEG (| log2FC | > 1) commonalities in our study were as high as 36.4% (Fig 3). We believe that this is meaningful because the commonalities among the previous studies [2,6,9,11] ranged from 6.3% to 60.6% (data were not shown). Additionally, commonalities with ASD were generally higher than those with other diseases (e.g., oral cancer or skin irritation), which would not expect an association with VPA exposure (S2(A) Fig) [43,44].
In our study, we found genes (adsl, mbd5, shank3a, and tsc1b) that had been reported to be ASD-related genes by multiple previous studies [2,6,9,11]. The adsl gene is involved in the de novo purine biosynthesis pathway [45]. Previous studies reported that deficiency of the adsl gene might result in infantile seizures and autism [45,46,47]. Although the exact mechanisms are not yet clear, the possible mechanisms might include deficient synthesis of purine nucleotides, impairment of the purine nucleotide cycle, and a build-up of defective enzyme substrates [45]. The mbd5 gene was thought to interact with myocyte enhancer factor 2c (MEF2C) in a complex bound to DNA [48]. Disruption of this complex due to mbd5 mutation altered gene expression in neurogenesis and neurodevelopment [48]. In a cohort study, deletion of mbd5 was also associated with autistic features [49]. The shank3 gene, which encodes synaptic scaffolding protein, is important for spine morphogenesis and synaptic plasticity [1,7,32]. It might play a critical role in the normal development of neuronal connectivity because its disruption at a genetic level causes defects at striatal synapses and cortico-striatal circuits [50], as well as development and speech delays [1,7], which might be associated with ASD. The tsc1 gene is thought to be linked to ASD via perturbation of cytoskeletal dynamics and dendritic spine structure by inhibition of the mTOR signalling cascade, which plays a crucial role in synapse protein synthesis [1]. Deletion of tsc1 in mice elicited a loss of cerebellar Purkinje cells and produced ASD-like behaviour, including abnormal social interactions and repetitive behaviours [51]. In our qPCR analysis, transcriptions of those genes were significantly altered (Fig 4). VPA, a chemical that induces ASD, led to a significant increase in the transcription of adsl5 and mbd5 but a significant decrease in the transcription of shank3 and tsc1b in our zebrafish embryo/larval model. These results were comparable with our results from RNA-seq analysis. As an exception, the direction of shank3a transcription was different. Among the exposed groups, however, the transcription pattern of shank3a was comparable. Despite the exception, the other results were considerably reliable.
Various genetic pathways, such as cell adhesion molecules, ion channels, scaffolding protein, cytoskeleton, and signalling pathways, are implicated in ASD [1]. We found some of the genes included in those pathways (nrxn2 (cell adhesion molecules), cacna1e, scn1a (ion channel), shank3 (scaffolding protein), actn4, gria3, grid1 (cytoskeleton), and eif4a (signalling pathway)) in our transcriptome analysis (S3-S7 Tables). Moreover, ASD-related genes such as ak1 (adenylate kinase 1) and ddb1 (damage-specific DNA binding protein 1) were included in the top 100 DEGs that we found in the 50 μM VPA-exposed group, based on the p-value (S10 and S11 Tables). Our DEG results imply that zebrafish might possess similar ASD genetic pathways to humans and other animal models. We found many of the common GO enrichments that are components of the ASD pathway [6,40]. Specifically, we reported 44 and 69 significantly altered GOs based on the transcription changes compared with the control group (Table 2) and concentration-dependent response (Table 3), respectively. The commonality was up to 40.0% and 29.4%, respectively, based on the p-values (p < 0.05 or p < 0.001). Those values were higher than the commonality with other diseases, i.e., oral cancer (8.7% and 1.0%, respectively) or dermatological disorder (19.2% and 1.9%, respectively) (S2(B) Fig) [52]. These comparison data of GO enrichments may support the idea that zebrafish, at least partly, share common physiology for ASD with other models.
We showed higher commonality with ASD than with other diseases, especially in the GO enrichment comparisons with stringent p-values (p < 0.001) (S2 Fig). However, commonality was relatively low in the DEG comparison. These results might be because GO enrichments represent a comprehensive approach. This implies that GO enrichments might adjust the bias from the transcriptions of certain genes. However, the results in S2 Fig are shown by a simple comparison among the respective lists of DEGs and GO enrichments, which could be due to methodological limitations. However, commonalities with other diseases, i.e., oral cancer or dermatological disorder, were observed in our study, although those were relatively weak. This might be expected because alterations of such basic pathways, such as the cytoskeleton or cell division, are also relevant to oral cancer or dermatological disorders [44,52]. In addition, VPA is known to have the potential to cause other side effects, such as hepatotoxicity and metabolic derangements, as well as VPA being an ASD inducer [53].
Among the significantly affected GOs, GOs involved in the cellular component pathway were abundant (Table 2). The presence of GO enrichments related to mitochondria (mitochondrion (GO: 0005739), mitochondrial inner membrane (GO: 0005743), and mitochondrial respiratory chain (GO: 0005746)) and the cytoskeleton (cytoskeleton (GO: 0005856) and actin cytoskeleton (GO: 0015629)) were significant (p < 0.001) ( Table 2). Recent studies have indicated that mitochondria dysfunction augments and disseminates brain abnormalities related to ASD [54,55]. Interestingly, a significant proportion of ASD patients have concomitant diseases related to mitochondrial disease and abnormal energy generation [56]. The cytoskeleton is also important for the positioning of cell adhesion molecules and neurotransmitter receptors at the synapse [2] and is essential for maintaining synaptic function [14]. The tsc1 gene, which we found to be significantly transcribed in our study, is related to the cytoskeleton pathway in ASD [1]. In addition, a previous study found neuronal signalling/cytoskeleton pathways containing actin cytoskeleton (GO: 0015629), actin cytoskeleton organization (GO: 0015629), and others to be major pathways of VPA-induced GOs in zebrafish when compared using the human functional network analysis [6] (S8 Table).
We found that significantly affected GOs were extended in a concentration-dependent manner (Table 3). We observed GO enrichments closer to the direct ASD pathway in VPAexposed zebrafish. GOs in the nervous system, e.g., synaptic transmission (GO: 0007268), transmission of nerve impulse (GO: 0019226), and neurological system process (GO: 0050877), were involved. Additionally, GO enrichments of synapses (synapse part (GO: 0044456) and synapse (GO: 0045202)) for which dysfunction is crucial for ASD were significant [2]. Furthermore, many of the ion channel-related GOs, which play an important role in the aetiology of ASD by enhancing neuronal excitability [1], were also significant, e.g., ion channel activity (GO: 0005216), cation channel activity (GO: 0005261), and ligand-gated channel activity (GO: 0022834). These significant changes observed in GOs in our study imply that zebrafish, at least partly, share commonalties with humans or mice with respect to ASD aetiology. This might indicate that zebrafish could be applicable to ASD research, e.g., ASD-inducing environmental chemical screening or ASD drug discovery.

Plausibility of zebrafish embryos/larvae as an ASD model based on behavioural changes
In addition to changes in gene transcription, the development of an experimental animal model for ASD requires strategies for behavioural observation. This is important because the symptoms of ASD are generally defined by behavioural characteristics, i.e., impaired social communication, repetitive behaviour, and cognitive deficits [25]. Previous studies have examined locomotor activity and social interaction behaviour in a transgenic monkey model [12]. Similar behavioural tests, such as social preference tests, the partition test, and the reciprocal social interaction test, have also been applied to rodent models [27,57]. Stewart et al. [25] compared measurable behavioural endpoints in rodent and zebrafish models and suggested a number of behavioural testing methods that could be applied to adult zebrafish. These included the social preference test, social interaction test, shoaling assessment, mirror stimulation test, and stereotypic circling and swimming [25]. However, adult zebrafish are not as effective as embryos for screening the ASD-like effects of many kinds of chemicals. The test tools for social behaviour in the zebrafish larvae stage have not yet been completely developed. Therefore, instead of social behaviour, locomotor activity and anxiety have been used as markers of the ASD phenotype in ASD toxicity tests involving zebrafish embryos [32,34,58]. Thus, in the Zebrafish embryos/larvae for autism present study, we analysed locomotor activity in zebrafish larvae after VPA exposure as a behavioural endpoint. We found that locomotor activity, i.e., the total distance moved and active duration, significantly decreased after exposure to 50 μM VPA (Fig 2). Previous studies reported similar levels of hypoactivity [31]. However, these effects, observed in both the current and previous studies, might be related to a VPA-induced delay in zebrafish hatching (Fig 1). We found that locomotor activity in zebrafish larvae slightly increased following exposure to 25 μM VPA. This was in accordance with a previous study and might indicate that VPA exposure causes anxiety in zebrafish larvae [34]. In a previous study, VPA exposure caused anxiety in larval zebrafish and resulted in social interaction deficits in adults [34]. Therefore, the observed effects on locomotor activity at the larval stage in our study could be linked to ASD-like behaviour. Although locomotor activity alone is not sufficient as a model of ASD in zebrafish larvae, this may be an important step, as such locomotor changes could reflect changes in basal activity elicited by ASD in terms of transcriptional alteration.

Further research regarding the use of zebrafish as an ASD experimental model
In the present study, we used transcriptome analysis and behaviour observation to show that the zebrafish embryo/larva has potential as an alternative animal model for ASD research. However, further studies are required. First, multiple chemicals used in ASD models, known as inducers, should be examined in the zebrafish embryo/larvae model. It is necessary to confirm that the zebrafish embryo/larvae model responds appropriately when other types of ASD inducers are used. Second, genes that are common to multiple ASD models and those that are particularly sensitive to environmental exposure should be determined. This genetic information may be useful in developing transgenic zebrafish models for screening the ASD-like effects of environmental chemicals and new drug candidates for ASD. Third, tools for appropriate behavioural examination should be optimized for the zebrafish embryos/larvae model, including, for instance, those that involve locomotor analysis.

Conclusion
In this study, we evaluated the plausibility of zebrafish embryos/larvae as an alternative animal model for ASD research. To this end, we exposed zebrafish embryos/larvae to VPA, examined the subsequent developmental and behavioural changes, and analysed transcriptional alterations by RNA-seq. We found behavioural and transcriptional changes that matched those from previous ASD studies using humans and rodent models. The results of this preliminary study indicate that the zebrafish embryo/larva has potential as an alternative experimental animal model for ASD toxicity tests and ASD drug discovery.