Zebrafish (AB line, standard length: 2.7–3.5 cm) were raised in a recirculating water system according to standard protocols (Westerfield, 2000). Aquaria-system water was dosed to a salinity of 500 μF with artificial ocean salt mix and buffered to pH 7.2 with NaHCO₃. After the group mating of four male and four female adult zebrafish (aged 6 months post-fertilization), embryos were collected and raised at 28.5°C in embryo medium (EM) at a density of 35–40 embryos per 10-cm diameter Petri dish. EM was prepared as previously reported (Coffin et al., 2009). Embryos were staged by hours post-fertilization (hpf) and days post-fertilization (dpf) as described previously (Kimmel et al., 1995). This study was carried out in accordance with the recommendations of the Fudan University Institutional Animal Care and Use Committee's guidelines (20120302-065). The protocol was approved by the Fudan University Institutional Animal Care and Use Committee.

Neomycin (5 g stock) or gentamicin (5 g stock; Sangon Biotech) was diluted in EM to final concentrations of 200 μM, or 25 and 50 μM, respectively. AG exposure paradigms were based on prior observations of the zebrafish lateral line (Owens et al., 2009; Stengel et al., 2017).

Healthy staged embryos were incubated from the 50% epiboly of embryonic shield stage (5.25 hpf) in neomycin or gentamicin to the long-pec stage during the hatching period (2 dpf). The medium was changed every 24 h. There was no statistical difference in embryo mortality between the exposed groups and control groups during this period. The 2 dpf zebrafish embryos were rinsed with EM several times and transferred to EM for further observation. The control groups were raised in EM the entire time, and the EM was changed every 24 h. Hair cell damage was examined in zebrafish embryos at 2 and 4 dpf.

Hair cells were pre-labeled with the mechanotransduction marker FM1-43 FX (3 μM in standard aquaria system water; Invitrogen Molecular Probes) for 45–60 s. The embryos were then quickly rinsed three times with EM. Using this procedure, FM1-43 FX is restricted to hair cells in neuromasts (Seiler and Nicolson, 1999). Hair cell survival was denoted by FM1-43 FX-positive cytoplasm surrounding the nucleus and an intact cell morphology.

Embryos were anesthetized with 0.001% MS-222 (Sigma). Zebrafish observations and image capture were conducted on a Zeiss Discovery V.20 microscope with a GFP filter set and an AxioCamHRc camera (Carl Zeiss), to visualize FM1-43 FX. Image stack projections or single image slices were exported from Slidebook software v. 4 (Olympus).

Total RNA for each sample was collected from 30 to 35 embryos at 2 dpf. Total RNA was isolated using Trizol reagent (Life Technologies). Agilent 2100 Bioanalyzer was used to characterize in vitro RNA transcripts for quality. The RNA integrity numbers (Rin) of all samples were higher than 8.0. The poly-A-containing mRNA molecules were purified using poly-T oligo-attached magnetic beads using two rounds of purification. A SuperScript Double-Stranded cDNA Synthesis kit (Invitrogen) was used to synthesize the double-stranded cDNA. Further library preparation was performed using TruSeq™ RNA Sample Preparation Kit (Illumina, cat# FC-122-1001). Libraries were sequenced as SR 2 × 100 bp using Illumina Hiseq2000 according to the manufacturer's instructions. We removed adaptor, low-complexity, and low-quality sequences in the raw reads. The remaining clean reads were used for further analyses.

Genomic DNA was isolated from the whole bodies of the four male and four female adult zebrafish from the mating pool using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's instructions. Equal amounts of total DNA from each zebrafish were pooled into a single sample. About 5 μg genomic DNA was sheared into fragments of 200–250 bp using a Covaris focused acoustic sonicator (Covaris, MA, USA). After size selection of fragments by 2% agarose gel electrophoresis, we constructed paired-end libraries with the NEXTflex DNA Sequencing Kit (BIOO Scientific, TX, USA) and Illumina adaptors using 1 μg of sheared input DNA. All libraries were sequenced on the Illumina Hiseq 2500 platform (Illumina Inc., CA, USA). After filtering out adaptor sequences, low-quality reads, and duplicate reads, a total of 54,802 Gb of data was retained for assembly.

The first 10 nts of paired-end reads were trimmed and filtered using a sequencing quality score cut-off of 30 for each RNA-seq dataset. Filtered reads were then aligned against the zebrafish reference genome (UCSC danRer7) with the Burrows-Wheeler Aligner program (BWA, version 0.7.12) using the default parameters (Li and Durbin, 2009). For each DNA-seq dataset, the first 10 nts of paired-end reads were also trimmed before being aligned against the reference genome with BWA. After mapping, the Samtools mpileup function was used to remove PCR duplicates, sort the alignment file, and call variants using default parameters for both RNA-seq and DNA-seq datasets (Li et al., 2009).

To improve the efficiency and accuracy of the discovery of RNA editing sites, we applied data quality filtration criteria during the three major steps of the analysis. Firstly, we took variant positions in the RNA into consideration if they varied from the reference genome and met our requirements in terms of the number, frequency, and quality of bases. We specifically required that the read coverage was at least 10, the mapping quality score of covered reads was at least 20, and there were at least two reads to support the alternative allele (Shamay-Ramot et al., 2015). Secondly, all homozygous and heterozygous DNA variants identified from DNA-seq data were filtered. Variants located within intergenic regions, polymers longer than 5 nts, or simple repeat regions, as well as known zebrafish SNPs in dbSNP (database version 135; http://www.ncbi.nlm.nih.gov/SNP/) and Ensembl variants annotated in Zv9 were also filtered. Thirdly, to further reduce false positives and identify RNA editing variants, we used the RNA-editing Calling process from a recently published method called SPRINT, which uses the editing type and cluster information of potential RNA editing sites to further enhance the calling power (Zhang et al., 2017).

To demonstrate differences in RNA editing among various treatments, all candidate RNA editing sites were combined as a candidate group, and then the sequencing statuses of these sites were re-evaluated with Samtools. For each different treatment, we calculated the ratio of edited sites for each replicated sample and determined the mean value of the treatment group. The degree of editing for a given site was calculated as the ratio of reads supporting the edited base to the total number of reads covering the site.

PCR amplification of gDNA and cDNA was carried out with the GeneAmp® PCR System 9700 (Thermo Fisher Scientific, MA, USA). The PCR protocol was as follows: 95°C for 4 min; 30 cycles of 94°C for 30 s, 57°C for 30 s, and 72°C for 30 s; 72°C for 10 min; followed by storage at 4°C as necessary. Direct sequencing was performed using an ABI 3730 DNA Analyzer (Applied Biosystems, MA, USA). The primers used for Sanger sequencing validation of called sites were listed in Supplementary Table 10. Sequence logos were generated with the program Two sample logo (Vacic et al., 2006).

Tophat v2.1.0 was used with default parameters to generate acceptable alignments for Cufflinks, in which the RNA-seq paired-end reads were aligned against the reference genome, UCSC danRer7 (Trapnell et al., 2012). Annotated gene expression was evaluated in FPKM (fragments per kilobase per million mapped fragments) using Cufflinks from RNA-seq data. The formula to calculate FPKM was as follows: FPKM = (number of mapping fragments) × 10³ × 10⁶/ [(length of transcript) × (number of total fragments)]. The method of normalizing expression data for comparison included log transformation and zero-mean normalization. Pathway analysis of genes with mRNA editing sites in zebrafish was performed using PANTHER v.10 (Mi et al., 2016).

The correlation between gene expression and editing was calculated using R software, version 3.1.2 (The R Foundation for Statistical Computing, Vienna, Austria; R Development Core Team, 2013).

For a given candidate RNA editing site, we used MIRANDA (Enright et al., 2003) to predict potential miRNA binding sites in both the genomic sequence and mRNA sequence after editing. Sequence windows of 50 bp around each editing site were provided to the prediction algorithm.

To establish a model of AG-induced ototoxicity, embryonic zebrafish were incubated in medium containing various concentrations of AG from the beginning of zebrafish ear development at the 50% epiboly/shield stage to 2 dpf, when the statoacoustic (VIIIth) ganglion becomes a separate section (Haddon and Lewis, 1996; Whitfield et al., 2002). Then, AG-treated zebrafish embryos were cleansed with and transferred to normal EM until 4 dpf, when the mechanoelectrical transducer channels attained maturity (Santos et al., 2006). Analysis of FM1-43 FX hair cell labeling at 2 and 4 dpf revealed that hair cells had differentiated into neuromasts in the head and along the posterior lateral line in the control zebrafish (Kimmel et al., 1995). However, few fluorescently labeled hair cell bundles were observed in the head and lateral line neuromasts of 2 and 4 dpf zebrafish embryos treated with AG (Figures 1A–H). The reduction in hair cells in the lateral line labeled with FM1-43X indicates the failure of hair cells to mature or to normally function under AG treatment.

To identify mRNA editing sites, we first sequenced mRNA from the 2 dpf embryos of the AB zebrafish line. We sequenced RNA from a total of nine AG treatment and control groups. Then, we sequenced DNA from adult zebrafish that represented the parental generation of the embryos. A total of 367.49 million mRNA sequencing reads and 219.2 million DNA sequencing reads were obtained. An average of 40.8 million reads or 4.0 Gb of mRNA-Seq data were generated per mRNA sample with a read length of 100 bp (paired-ends) and expected insertion size of 200 bp. The alignment rate against the reference genome (see section Materials and Methods) ranged from 94.04% (sample Genta_50_rep2) to 95.78% (sample Genta_50_rep1), with average mapping rates of 95.06 and 94.83% for control and drug-treated samples, respectively (Supplementary Table 1).

We identified an average of 6,973 edited RNA sites in gene body regions per RNA sample, with the number of edited sites ranging from 6,859 to 7,019. Exposure to 25 μM gentamycin resulted in hair cell loss similar to that of 50 μM gentamycin; therefore, we combined the data obtained from these samples in subsequent analyses. Finally, we identified 6,850 overlapping mRNA editing sites and calculated the level of editing in control, gentamycin, and neomycin-treated samples (Supplementary Table 2). Of the edited sites, 5,662 consisted of canonical, A-to-G or T-to-C editing. Two types of transitions, G-to-A and C-to-T, accounted for 77.4% of the non-canonical events observed.

We verified a subset of the edited RNA sites using PCR amplification and Sanger sequencing of DNA and mRNA (reverse-transcribed to cDNA). We examined 182 editing sites across 17 genes, including 13 G-to-A and C-to-T editing events. The validation results yielded a false-discovery rate of 5.33% (9/169) for A-to-G or T-to-C sites identified in our analysis, which was lower than the false-discovery rate for canonical sites in the previous study (Peng et al., 2012). The false-discovery rate for the non-canonical type was ~84.6% (11/13; Supplementary Table 3). The validation of editing events in the sp4 and si:ch211-114n24.7 genes is shown as representative examples (Supplementary Image 1).

We identified 6,850 mRNA editing sites distributed over 762 genes. About 80% of these sites were A-to-G substitutions. The master list of editing sites contains 5,590 exonic sites, 1,243 intronic sites, and 17 sites located in regions with conflicting database annotations. Among exonic sites, the UTRs, especially the 3′ UTR, contained the greatest percentage of both A-to-G and non-A-to-G variants. CDS regions contained a total of 479 edited sites, of which 35% led to amino acid changes. Furthermore, there were significantly more non-A-to-G sites than A-to-G sites in CDS regions (p < 2.2 × 10⁻¹⁶). Among the 762 genes with editing sites identified in this study, 299 genes were reported to contain editing sites in previous studies of the human transcriptome (Peng et al., 2012; Ramaswami et al., 2013; Qiu et al., 2016), 54 genes had been described as edited in the mouse transcriptome (Danecek et al., 2012; Gu et al., 2012), and 325 genes had been described as edited in the zebrafish transcriptome (Shamay-Ramot et al., 2015). There were 279 genes which had not been described as edited before.

We also found that the RNA editing sites were clustered, with an average of nine editing sites per gene. The extent of A-to-G site clustering in our data set, with 60.06% of sites arranged in clusters of ≥3 sites within 100 bp, is lower than what is found in the DARNED database (85.02%) but is higher than that observed in another deep-sequencing data set acquired by Peng et al. (30.89%) (Kiran and Baranov, 2010; Peng et al., 2012). A large number of transcripts presented multiple editing sites (204 genes with ≥10 sites each), such as sp4, slc22a31, plekhm1, adrbk2, pkn1, and samhd1, which were validated with PCR and Sanger sequencing.

Pathway analysis of genes with mRNA editing sites in zebrafish was performed using PANTHER v.10 (Mi et al., 2016). mRNA editing sites were significantly enriched (p < 0.05) in the Huntington disease, semaphorin-mediated axon guidance, angiogenesis, pyrimidine metabolism, cytoskeletal regulation by Rho GTPase, and Alzheimer disease-presenilin and FGF signaling pathways (Supplementary Table 4). Moreover, plch1 and snap29, part of the 5HT2-type receptor-mediated signaling pathway, also contained several editing sites. These findings indicate that RNA editing participates in the maintenance and regulation of zebrafish development, and particularly in neural development.

RNA editing can affect anywhere from 0 to 100% of an RNA population and leads to a mixed RNA population. To explore the influence of AG on the mRNA editome, we investigated differences in the degree of editing between AG-treated and control samples (Figure 3A). Editing sites with more than a 15% change in the degree of editing in both AG treatments were included in the high variation group and used for further analysis. Notably, in the high variation group, 687 sites spanning 333 genes accounted for 43% of all detectable RNA-edited genes (Supplementary Table 5). Of these RNA editing sites, 462 and 225 showed a lower and higher degree of editing, respectively, in AG-treated samples than in control samples. In the high variation group, canonical variants (A-to-G or T-to-C) accounted for 75.8% of all RNA edits. The average value of the editing degree variation was about 26% among down-regulated edits and 27% among up-regulated edits. A total of 21 editing events across 16 genes were likely to change the encoded amino acids, and 57% of these changes were non- canonical variants. Pathway analysis showed that the 333 genes in the high variation group were significantly enriched (p < 0.05) in axon guidance mediated by netrin and in the FGF and Ras signaling pathways (Supplementary Table 6).

A total of 178 genes with editing events were significantly enriched in lateral line cells or lateral line placode cells, with 76 genes containing editing events with editing variations of more than 15% between control and AG-treated samples (Jiang et al., 2014; Steiner et al., 2014; Figure 2). Moreover, 24 editing sites were identified as having more than a 15% change in the level of editing following AG treatment in 12 genes potentially related to ototoxicity, including plekhm1, fgfr1a, sox9a, and calrl2 (Karasawa et al., 2011; Girotto et al., 2013; Stamatiou and Stankovic, 2013; Azaiez et al., 2014; Table 1). Editing events in the deafness-related genes bdp1 and alms1 led to amino acid changes, which may be one of the potential pathways of AG-induced ototoxicity. Other important genes with editing events with more than 15% variation that may change the amino acid sequence are listed in Table 2. Furthermore, editing variation in the mitochondrial inner membrane-expressed gene oxa1l could also affect the biological function of the protein.

We wanted to investigate whether there were sequence characteristics underlying the AG targeting of RNA. Therefore, we compared the flanking sequences of editing sites with editing variations >15% of control and AG-treated samples (high variation group) with those that had <5% variation in the degree of editing (low variation group; Figure 3B). We found that sites in the high variation group were flanked by G-enriched sequences and exhibited a depletion of C nucleotides downstream of the edited site. This agrees with previous findings that the C nucleotide is underrepresented at the +1 position of editing sites in zebrafish, while G is overrepresented at the +1 position of editing sites in humans and mice (Shamay-Ramot et al., 2015). This is also in agreement with the 3′ nearest neighbor preferences of the ADAR2 enzymes, with G being preferred (Eggington et al., 2011). Additionally, the presence of a T nucleotide four bases downstream of the edited site was characteristic of sites with down-regulated degrees of editing (more than 15%) in AG-treated samples. Enrichment of G nucleotides three bases upstream of the edited site was characteristics of sites with up-regulated degrees of editing (more than 15%) in AG-treated samples (Figure 3C). Nucleotide preferences in the flanking sequences of RNA editing sites indicate that the binding of AG is facilitated by specific sequence and structural features.

We measured variations in gene expression between the AG-treated and control 2 dpf embryo samples. RNA of all samples were isolated from whole embryos. The expression of 988 genes differed more than 1.5-fold between both neomycin and gentamycin-treated samples and normal control samples. Of these genes, 257 genes were lateral line-enriched genes identified in previous studies (Jiang et al., 2014; Steiner et al., 2014). Among the 722 genes with down-regulated expression in the two AG-treated samples (Supplementary Table 7), 113 genes were also down-regulated more than 1.5-fold in the expression analysis of lateral line cells at 1 h after neomycin-induced hair cell death, including her6, etv4, and tcf7l1a (Jiang et al., 2014). Twenty-seven of the down-regulated genes, including ptprq, ush1g, climp-63, and genes involved in the Wnt, Notch, and FGF signaling pathways, play important roles in ear development or ototoxicity (Karasawa et al., 2010; Stamatiou and Stankovic, 2013; Jiang et al., 2014; Table 3). There were 266 genes with up-regulated expression in AG-treated samples, among which 47 genes were up-regulated more than 1.5-fold in the lateral line cells at 1 h after neomycin-induced hair cell death (Jiang et al., 2014). Pathway analysis showed that down-regulated genes were significantly enriched (p < 0.05) in cytoskeletal regulation by Rho GTPase and the thyrotropin-releasing hormone receptor, cadherin, and Wnt signaling pathways (Supplementary Table 8). Up-regulated genes were significantly enriched (p < 0.05) in DNA replication and in the apoptosis signaling pathway (Supplementary Table 9).

The RNA-seq data allowed us to assess the connection between RNA editing and gene expression in the experimental samples. Increases in the level of RNA editing were weakly, but significantly, correlated with increases in gene expression in the gentamycin treatment group (Pearson's correlation test: r = 0.10, 95% CI: 0.06–0.20, p = 0.0088; Figure 4). Moreover, 160 editing sites were located in 19 genes with more than a 1.5-fold difference in expression between AG-treated and control samples. The degree of editing changed more than 15% in 21 of these editing sites, and all but two were located in non-coding regions (Table 4). Among these 21 editing sites, three were predicted by MIRANDA to alter or create miRNA binding sites and possibly affect the expression of shmt2, si:ch211-241b2.1, and plekhm1 (Table 5). Under AG treatment, the expression levels of shmt2 and plekhm1 were down-regulated more than 1.5-fold with the up-regulated editing of these two miRNA binding sites. shmt2 (serine hydroxymethyltransferase 2) is mainly distributed in the mitochondria and plays an important role in mitochondrial DNA synthesis and glycine production, while plekhm1 is associated with syndromic hearing loss (Van Wesenbeeck et al., 2007; Anderson and Stover, 2009).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.