Transcriptome profiling of zebrafish optic fissure fusion

Incomplete fusion of the optic fissure leads to ocular coloboma, a congenital eye defect that affects up to 7.5 per 10,000 births and accounts for up to 10 percent of childhood blindness. The molecular and cellular mechanisms that facilitate optic fissure fusion remain elusive. We have profiled global gene expression during optic fissure morphogenesis by transcriptome analysis of tissue dissected from the margins of the zebrafish optic fissure and the opposing dorsal retina before (32 hours post fertilisation, hpf), during (48 hpf) and after (56 hpf) optic fissure fusion. Differential expression analysis between optic fissure and dorsal retinal tissue resulted in the detection of several known and novel developmental genes. The expression of selected genes was validated by qRT-PCR analysis and localisation investigated using in situ hybridisation. We discuss significantly overrepresented functional ontology categories in the context of optic fissure morphogenesis and highlight interesting transcripts from hierarchical clustering for subsequent analysis. We have identified netrin1a (ntn1a) as highly differentially expressed across optic fissure fusion, with a resultant ocular coloboma phenotype following morpholino antisense translation-blocking knockdown and downstream disruption of atoh7 expression. To support the identification of candidate genes in human studies, we have generated an online open-access resource for fast and simple quantitative querying of the gene expression data. Our study represents the first comprehensive analysis of the zebrafish optic fissure transcriptome and provides a valuable resource to facilitate our understanding of the complex aetiology of ocular coloboma.


RNA collection and sequencing
Eye diameters at 32, 48 and 56 hpf were estimated by imaging fifty staged embryos from a single clutch. Horizontal eye diameter was measured in the developing eye using ImageJ. Ten such clutches per timepoint were collected and measured independently and the population means of the eye diameters, estimated from the means of clutch replicates, were determined to be 186.1 μm for 32 hpf (95% CI = 178.0 μm-194.1 μm), 212.7 μm for 48 hpf (95% CI = 206.7 μm-218.6 μm) and 232.4 μm for 56 hpf (95% CI = 223.5 μm-241.3 μm). During sample dissection, eye diameter was measured and only those that fell within the 95% confidence intervals of the population mean, and fulfilled conventional staging criteria, were chosen for dissection. Staged embryos were fixed for 1-2 hours by immersion in RNase-free 4% paraformaldehyde at 4°C. Optic fissure (OF; ~50 µm from either side of the fissure) and opposing dorsal retinal tissue (DR; ~100 µm) was dissected using fine tip tweezers and a sapphire blade (World precision instruments), snap-frozen and stored for RNA extraction.
Tissue was dissected from 5 biological replicates per timepoint; each biological replicate was generated from an independent breeding pair. RNA was extracted using the RNeasy FFPE Kit (Qiagen) and quantified using the Bioanalyzer 2100 RNA Pico system (Agilent biosystems) per manufacturer's instructions. Samples of RIN ≥ 8 were selected for sequencing. cDNA libraries were constructed from total RNA using the Clontech SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Clontech) and sequenced on a HiSeq 2500 system using v4 chemistry (Illumina). Libraries were paired-end sequenced to a depth of between 20-40 million reads at a length of 100 bp. Library preparation and sequencing were performed by Eurofins Genomics (Ebersberg, Germany).

RNA-seq data analysis
Raw sequencing reads per sample per lane were concatenated and read quality assessed.
Initially samples were evaluated with FASTQC (v0.11.5) demonstrating between 40-80% duplication of reads per sample. Samples with extreme duplication rates were removed from the analysis. Adapter and low quality scored base reads (Phred < 6) were removed using Trimgalore (v0.4.3). High-quality cleaned FASTQ reads were aligned to the Ensembl D.rerio reference genome (build GRCz10) and annotation (build GRCz10.89) downloaded June 2017 using the splice-aware aligner STAR (v2.5.2b). Resultant BAM files were sorted, indexed, duplicates marked with Picard tools and subsequently used to generate count files using HTSeq (v0.61, 1 ). To assess mapping efficiencies, aligned reads were analysed with RNA-SeQC (v1.1.8, 2 ). Statistics of RNA-seq reads and mapping were outlined in supplementary file SI_01. Duplicate reads were not removed from the count data, however D.rerio rRNA and tRNA species were. Count files generated using HTSeq were processed for subsequent differential expression analysis with the Bioconductor R packages tximport 3 and DESeq2 (v1.18.1, 4 ). PCA performed on the regularised log transformed count data in DESeq2.
To investigate differential gene expression we performed successive pairwise condition testing (i.e. OF tissue 32 hpf vs OF tissue 48 hpf, OF tissue 48 hpf vs OF tissue 56 hpf, and so on) using the raw count data. DESeq2 was used to model the counts to the negative binomial GLM, normalising to library size (factor), and identify differentially expressed gene Visualisation of DEG events was carried out using the Bioconductor package NMF 6 .
Significant DE Ensembl gene ID lists were generated and heatmaps produced through clustering using hierarchical average linkage clustering and Euclidean distances.

Gene Ontology Analysis
Gene Ontology (GO) analysis of DEGs was performed using a combination of Bioconductor packages; GOseq 7 , GenomicFeatures 8 and Pathview 9 . GOseq was used with correction for gene length bias in the RNA-seq data. The GO hierarchy "Biological process" was used to query these genes against the background of all genes expressed the RNA-seq dataset. Only GO terms with an adjusted p < 0.05 were considered. Redundant GO terms were removed with REViGO 10 .

Quantitative real-time PCR analysis (qRT-PCR)
For verification of transcript changes identified in RNA-seq, OF and opposing DR tissue was dissected from ten 32, 48 and 56 hpf embryos (one eye per embryo) and pooled for RNA extraction. Each pool was defined as one of three biological replicates, obtained from independent breeding pairs. RNA was extracted using the RNeasy FFPE Kit (Qiagen), quantified using the RNA Analysis 2200 ScreenTape system (Agilent Genomics) and Co-injection of p53 morpholino was carried out with the gene-specific morpholinos as control experiments.