Transcriptome analysis of the zebrafish atoh7-/- Mutant, lakritz, highlights Atoh7-dependent genetic networks with potential implications for human eye diseases.

Expression of the bHLH transcription protein Atoh7 is a crucial factor conferring competence to retinal progenitor cells for the development of retinal ganglion cells. Several studies have emerged establishing ATOH7 as a retinal disease gene. Remarkably, such studies uncovered ATOH7 variants associated with global eye defects including optic nerve hypoplasia, microphthalmia, retinal vascular disorders, and glaucoma. The complex genetic networks and cellular decisions arising downstream of atoh7 expression, and how their dysregulation cause development of such disease traits remains unknown. To begin to understand such Atoh7-dependent events in vivo, we performed transcriptome analysis of wild-type and atoh7 mutant (lakritz) zebrafish embryos at the onset of retinal ganglion cell differentiation. We investigated in silico interplays of atoh7 and other disease-related genes and pathways. By network reconstruction analysis of differentially expressed genes, we identified gene clusters enriched in retinal development, cell cycle, chromatin remodeling, stress response, and Wnt pathways. By weighted gene coexpression network, we identified coexpression modules affected by the mutation and enriched in retina development genes tightly connected to atoh7. We established the groundwork whereby Atoh7-linked cellular and molecular processes can be investigated in the dynamic multi-tissue environment of the developing normal and diseased vertebrate eye.

eye field transcription factors

retinal non‐attachment

optic nerve aplasia

optic nerve hypoplasia

persistent hyperplastic primary vitreous

retinal ganglion cells

retinal progenitor cells

INTRODUCTION

Retinal ganglion cells (RGCs) collect visual information from the neural retina in the eye and convey it to the visual cortex of the brain. In healthy people, this information is transmitted along the optic nerve, which is composed mainly of axons formed from the cell bodies of RGCs. Inherited diseases affecting the development of RGCs and the optic nerve can interrupt this information flow causing permanent blindness. , Atoh7 is an evolutionarily conserved, developmentally regulated transcription factor crucial for the genesis of RGCs in different vertebrate models. , , , , Studies have shown that induced or naturally occurring mutations in the atoh7 gene result in retinal progenitor cells (RPCs) failing to develop into RGCs and the optic nerve. , , Likewise, an increasing number of studies highlight ATOH7 as an emerging candidate for eye diseases in humans. Variations in the ATOH7 locus have been associated with optic nerve hypoplasia (ONH) and aplasia (ONA), , , , further pointing towards the crucial role of atoh7 in RGC genesis and optic nerve development. Remarkably, a number of studies have also emerged, which highlight ATOH7 variants as associated with multiple eye disease traits. These include disorders of the retinal vasculature, such as retinal non‐attachment (NCRNA) and persistent hyperplastic primary vitreous (PHPV) (OMIM:# 221900, ORPHA:91495), , , , , , , , as well as corneal opacity, microcornea, microphthalmia (ORPHA:289499), , and glaucoma. The Atoh7‐regulated gene networks involved, and how their disruption contribute to the development of such ocular disorders remain unknown.

The zebrafish has long been valued as a paradigm for disentangling the genetics and cell biology of fundamental eye developmental processes. , The rapidly and externally developing transparent zebrafish embryos are amenable to easy genetic manipulation, allowing fast generation and identification of mutants modelling human ocular genetic disorders. , , , , , , Such disease models can be concurrently investigated in large‐scale genetics, drug screening, in vivo cell biology of early disease development as well as behavioral assays. , , , These potentials substantially aid fast progress in the validation of human genome association studies and in preclinical therapy development paths toward the early diagnosis and/or restoration of visual function. , , ,

We here begin to explore the potentials of the lakritz zebrafish mutant carrying a loss of function mutation in the atoh7 gene. With the analysis of available microarray data, we investigate Atoh7‐regulated gene networks and interrogate how deregulation of these networks during early onset of RGC genesis might contribute to the development of atoh7‐associated eye disorders. We provide a cohort of statistically significantly regulated Atoh7 target genes, including previously known Atoh7 targets such as atoh7 itself. , Remarkably, at this early RGC developmental time‐point, the most significant targets comprehend previously unreported eye field transcription factors, Wnt signaling pathway components, chromatin and cytoskeletal regulators, and even stress‐response proteins as major Atoh7‐regulated genes. Furthermore, several components of these pathways include eye disease gene markers.

With these data in hand, we can now begin to exploit the power of zebrafish as in vivo vertebrate model to assess how dysregulation of one or more of these components might affect the coordinated development of ocular tissues. This will hopefully provide us with a next step forward in the identification of potential targets for the early detection and/or specific treatment of inherited eye diseases such as retinal‐vascular disorders.

MATERIALS AND METHODS

Wild‐type and transgenic zebrafish

Fish used in this study were identified heterozygous carriers of the lakritz mutation crossed in the (AB/AB) background as well as transgenic tg(lakritz/atoh7:gap43‐RFP) heterozygous carriers. , All fish were maintained at 26°C‐28°C as described previously. , Embryos were obtained by breeding adult male and female fish at ratio 1:1. After fertilization, eggs were collected and maintained at 28.5°C and staged using standard morphological criteria. Fish were kept and experiments were performed in accordance to local animal welfare agencies and European Union animal welfare guidelines.

Eyes and body sample collection

Single pairs of eyes were dissected from single embryos at 25, 28, 35, 48, 72, and 96 hpf and snap‐frozen in liquid nitrogen, and stored at −80°C. Embryos older than 28 hpf were first anesthetized for 5‐10 minutes in ethyl 3‐aminobenzoate methanesulfonate (MS‐222) (Sigma‐Aldrich) in E3 medium. The corresponding body of each embryo was collected and used immediately to perform genotyping analysis to identify the corresponding lakritz and wild‐type eyes. All embryos were collected from the same batches of fish stock to maintain a uniform genetic background.

DNA extraction from zebrafish body biopsies and genotyping

Genomic DNA extraction from each single body was performed in 100 μL of lysis buffer containing Proteinase K‐20 mg/mL (EuroClone S.p.A. Milan, Italy), 2 mol/L Tris‐HCl ph 8.0, 0.5 mol/L EDTA ph 8.0 and 5 mol/L NaCl, 20% SDS in a final volume of 50 μL ultra H20. After 3 hours of incubation at 65°C, the gDNA was purified with an ethanol precipitation step and re‐suspended in 50 μL of Dnase/Rnase H2O. The genotyping was performed by restriction fragment length polymorphism (RFLP) assay as previously described. An ∼300 bp fragment of atoh7 was PCR amplified with 1 U Taq DNA polymerase (Applied Biosystems by Life Technologies, Carlsbad, CA, USA) according to the manufacturer protocols in a 30 μL PCR mix containing 100 ng of purified gDNA (from each single embryo body) with the following primers: Forward 5′‐CCGGAATTACATCCCAAGAAC‐3′ and Reverse 5′‐GGCCATGATGTAGCTCAGAG‐3′. The PCR conditions were as follows: initial denaturation (95°C for 5 minutes), followed by 40 cycles of denaturation (95°C for 45 seconds), annealing (56°C for 45 seconds), extension (72°C for 45 second), and a final extension at 72°C for 5 minutes. The resulting PCR product was digested with StuI restriction enzyme (NEB), according to the manufacturer's protocols. The digested product was analyzed on a 2% agarose gel in 1X Tris‐Acetate EDTA (TAE) buffer (Sigma‐Aldrich) to highlight wild type or lakritz mutant corresponding fragments. The L44P mutation eliminates a restriction site found in the L44 allele and therefore can be visualized as an undigested ∼300 bp fragment rather than the ∼100 and ∼200 bp fragments expected from the wild‐type condition (Figure 1B). 1Kb DNA Ladder was used as a reference (Gene ruler 1Kb plus, Thermo Scientific).

FIGURE 1

Scheme of the experimental design for the comparative array analysis. A, Confocal images showing examples of wild‐type and lak−/− (lakritz); tg(atoh7:gap43‐RFP) embryos at 96 hpf. The RFP‐positive optic chiasm and RGCs are absent in the retina of a lakritz embryo. B, Pairs of eyes were dissected from single embryos at 25, 28, 35, 48, 72, and 96 hpf . Genotyping on the gDNA extracted from each corresponding cell body was performed to identify lakritz and wild‐type embryos (see materials and methods section). The RNA extracted from each pair of eyes corresponding to either a lakritz or wild‐type embryo was amplified and used for the microarray analysis and qRT‐PCR expression analysis. [Color figure can be viewed at wileyonlinelibrary.com]

Affymetrix arrays hybridization and analysis

For the microarray analysis, three pairs of wild‐type and lakritz 28‐30 hpf embryos representing three biological replicates were dissected and placed in Trizol reagent (Thermo scientific Life Technologies) for the total RNA extraction according to the manufacturer's instructions. T7‐based linear amplification of the mRNA was performed using the megascript kit from Ambion. Hybridization was performed on the Affymetrix GeneChip platform and processed according to standard procedure.

Correspondence between Zebrafish Affymetrix probesets and EnsEMBL gene annotations was retrieved using BioMart (EnsEMBL Version 84, March 2016). Batch effect removal was applied to adjust for known batch effect by first filtering the normalized matrix of intensities discarding probes with total abundance between samples lower than first quartile (Q1 = 16.26), and then correcting intensities using ComBat package with extraction day as the known batch (Figure S1). Differential expression analysis between mutant and wild‐type samples was performed using Limma package. For probes showing statistically significant differential expression (adj. P value < 0.05), annotations of corresponding genes were retrieved from the EnsEMBL database using BiomaRt package.

Quantitative Real‐Time PCR (qRT‐PCR)

For the qRT‐PCR analysis of anillin and atoh7, total RNA from five pulled pairs of frozen eyes, corresponding to either lakritz or wildtype embryos, was used. After Turbo DNase treatment (Thermo scientific‐Ambion), according to the manufacturer instructions, the RNA concentrations were measured with a Nanodrop ND‐1000 spectrophotometer (NanoDrop Technologies Inc). The RNA integrity was verified by loading the samples on 1% agarose gel that was run at 100 Volts in TBE 1X. 500 ng of extracted RNA from each sample was retrotranscribed with a RevertAidTM First Strand cDNA Synthesis Kit (Thermo SCIENTIFIC), following the manufacturer's protocol. Quantitative real‐time PCR reactions were performed on a Bio‐Rad CFX96 Termo‐cycler with Kapa Syber Fast qPCR master mix (2×) kit (Sigma‐Aldrich), according to the manufacturer's instructions. Templates were 1:10 diluted cDNA samples. For the negative controls, cDNAs were replaced by DEPC water. All real‐time assays were carried out using 10 ng of cDNA. The PCR profile was as follows: 15 seconds at 95°C, followed by 40 cycles 60°C for 20 seconds, 72°C for 40 seconds. For the melting curve, 0.5°C was increased every 5 seconds from 65°C to 95°C. All reactions were run in triplicate and both glyceraldehyde‐3‐phosphate dehydrogenase (GAPDH) and Ubiquitin Conjugating Enzyme E2A (UBE2A) were used as reference genes. Each experiment was performed in triplicate and repeated two times. The relative gene expression was calculated using the ΔΔCT method. Statistical analyses were performed with Prism 5 (GraphPad Software), and statistical significance was set to P < .05 for all experiments. The values are expressed as mean ± SEM, and the differences between groups were investigated using unpaired two‐tailed Student's t test (GraphPad Software). A list of primers is shown in Table 1.

TABLE 1

List of primers with amplicon sizes used for quantitative real‐time PCR

Primers	Sequence (5′> 3′)	Product length (bp)
GAPHD_zf FOR	TCACAAACGAGGACACAACCA	219
GAPHD_zf REV	CGCCTTCTGCCTTAACCTCA
Ube2a_zf FOR	CTGAAGGAACACCTTTTGAAGATG	215
Ube2a_zf REV	GATCCAGTAAAGACTGTATTGAG
Atoh7_zf FOR	TCACCTGTGGAAAGTGACTG	254
Atoh7_zf REV	CTCATTCACAACCCGCCCAA
Anln_zf FOR	AAAGGCTTCCTGACTATGTTTG	107
Anln_zf REV	CATCATCAGGGTAGGTCCA

Functional category enrichment analysis and network analysis

Functional enrichment analysis was performed with Metascape using the Danio rerio EnsEMBL IDs of the list of differentially regulated genes as “Input as species” and “Analysis as species” species through the custom analysis mode. Enrichment analysis was performed against GO “Biological Process” using P value cut‐off of < 0.05 and otherwise default parameters. Human disease annotation was performed with Metascape using Danio rerio as “Input as species” and H sapiens as “Analysis as species.” Under the “Annotation” mode, all repositories under Genotype/Phenotype/Disease were selected for disease annotation. To rule out potential sampling or biological bias, we performed an additional enrichment analysis by restricting the background to the list of expressed genes as detected by the arrays, using both Metascape and KOBAS. Network interactions between the differentially expressed genes were retrieved through the STRING database, “multiple proteins” mode and default parameters otherwise. Network interactions were visualized using Cytoscape.

Weighted gene co‐expression network analysis

The pipeline proposed by Langfelder and collaborators in their CRAN package was followed to infer gene co‐expression networks and identify network modules within R 3.6.3 statistical environment. Networks were inferred using the TOMsimilarityFromExpr function with “cor” as gene coexpression measure. The soft‐threshold parameter was optimized with the function pickSoftThreshold and the best threshold (α = 16) selected by visual inspection to follow a scale‐free topology model, as suggested by the WGCNA pipeline. Correlations between modules eigengenes, status, and library batch were computed. Modules with the highest correlation for status and significant P value (α ≤ 0.05) were selected for further analysis. Within the selected modules, highly connected structure of submodules were identified using the leading eigenvector community detection method implemented in the igraph package (v1.2.5) for R (http://igraph.com).

Whole mount immunohistochemistry

For immunohistochemical labelling, embryos were fixed in 4% PFA for 1 hour at room temperature or overnight at 4°C. Embryos were washed three times in PTw and kept for a week maximum in PTw. Fixed embryos were blocked in blocking solution (10% goat serum, 1% bovine serum albumin, and 0.2% Triton X‐100 in PBS) for 1 hour. Embryos were permeabilized with 0.25% trypsin‐EDTA (1X, Phenol red, Gibco; Life Technologies) on ice for 5 minutes. Primary (mouse anti‐β‐catenin, 1:100, Cat. No 610153, BD Biosciences) and secondary (anti‐mouse IgG conjugated to Alexa Fluor 488, 1:250, Life Technologies) antibodies were added for two overnights each and DAPI was added from the first day of incubation in the antibody mix. Stained embryos were kept in PTw at 4°C in dark until imaging. Embryos were embedded onto a 35 mm Glass‐bottom Microwell dish (p35G‐1.5‐10‐C, MatTek) and oriented with a femtoloader tip (Eppendorf, Leipzig, Germany) in the position needed for imaging until the agarose had polymerized. Confocal imaging was performed using a laser scanning confocal microscope Leica SpE using a Leica 40×, 1.15 NA oil‐immersion objective.

For the analysis of the intensity along the apical‐to‐basal membrane, a line of a defined length along the apical‐to‐basal membrane was drawn for nine cells at three different confocal z‐sections for each embryo. Signal intensities were obtained using Fiji and average values for the nine cells were calculated. For the analysis of the intensity along the apical membrane, the number of peaks were counted after Ctnnb1 signal intensity measurement. For the signal intensity measurements, a line of defined length was drawn along the apical membrane of the retina on Fiji and signal intensities were retrieved. The length of line was the same for all z‐sections of an individual embryo. Normalization was performed by the highest value for each line. Measurements were obtained on three different Z sections for each embryo.

RESULTS

Transcriptome analysis of wild‐type and lakritz identified 137 statistically significant differentially expressed genes

Expression of atoh7 in the retina is first detected at around 25‐28 hpf and it reaches its peak at around 36 hpf. The earliest post‐mitotic RGCs in the retina are detected at around 28 hpf, a developmental time‐point corresponding to the earliest onset of retinal differentiation. To identify early Atoh7‐regulated genes, transcriptome analysis was performed on eyes from single lak mutant (lakritz) and wild‐type zebrafish embryos at 28‐30 hours post‐fertilization (hpf) based on Affymetrix microarrays (see methods and Figure 1).

Differential gene expression analysis performed with LIMMA resulted, after batch effect removal, with 171 significantly differentially expressed probes (adj. P value < .05) (Figure 2A) corresponding to 137 genes annotated onto EnsEMBL database (Figure 2B and Table S1). Among these, we confirmed downregulation of atoh7 in the lakritz, consistent with its known role as self‐activator. Also consistent with the presence of bona fide Atoh7‐regulated targets in the 137 cohort is the presence of additional seven genes, which have been previously reported to contain a well‐characterized Ath5 consensus binding site. , Among the eight Atoh7‐direct targets, besides atoh7 itself, the thyrotroph embryonic factor tefa and the atypical cadherin receptor 1 celsr1a (CELSR1) were downregulated in the lakritz, suggesting their positive regulation by Atoh7. Conversely, the retina and anterior neural fold homeobox transcription factor rx1 (RAX), the Wnt signaling pathway regulator notum, the transmembrane protein tmem165, and the F‐actin binding protein and cytokinesis regulator anillin (ANLN) were upregulated in the lakritz, suggesting that they are negatively modulated by Atoh7.

FIGURE 2

Volcano and heatmap of differentially expressed genes in lakritz vs wild‐type eyes A,Volcano plot highlighting Atoh7 and its direct targets (in red) among all differentially expressed probes (in green) with adjusted P value < 0.05. B, Heatmap was constructed by calculating row Z‐score using normalized log2 intensities of 144 of the 171 differentially expressed probes with corresponding gene annotation, using complete hierarchical clustering in R. [Color figure can be viewed at wileyonlinelibrary.com]

Functional category enrichment reveals neural retina, cell cycle and Wnt pathways regulated downstream of Atoh7

Functional enrichment analysis in Metascape using default parameters reveals “neural retina development” (GO:0003407) as the most highly significantly enriched GO Biological Process category, consistent with the role of Atoh7 as regulator of retinal development (Figure 3 and Table S2). To further increase stringency, ruling out sampling or biological bias in the analysis, the background was restricted to the list of expressed genes as detected by the arrays using both Metascape and Kobas 3.0, which incorporates knowledge from five pathway databases (KEGG PATHWAY, PID, BioCyc, Reactome, and Panther) and five human disease databases (OMIM, KEGG DISEASE, FunDO, GAD, and NHGRI GWAS Catalog). This analysis consistently underscored “neural retina development” (GO:0 003 407) as the enriched term both by Metascape (not shown) and Kobas (Table S3). This category contains nine statistically significantly differentially expressed genes (including atoh7), which comprehend early expressed eye‐field transcription factors (EFTFs, rx1, and six6a), , , , stress response and extracellular matrix remodeling factors (hsp70.1, mmp14a), , , , chromatin regulators (smarca5) as well as microtubules organizers and cell cycle regulatory proteins (tubgcp4, znf503, and gnl2) , , (Table 2). Other (albeit less significant) relevant biological processes emerging from this analysis were “cell cycle process” (GO:0022402), “chromatin remodeling” (GO:0006338), and “Wnt signaling pathway” (GO:0016055) (Figure 3 and Table S2).

FIGURE 3

Functional enrichment analysis. Statistically significantly over‐represented GO Biological Process categories (Metascape). See also Table S2 [Color figure can be viewed at wileyonlinelibrary.com]

TABLE 2

Significantly differentially expressed genes belonging to the “neural retina development” category (see also Figure 3 and Table S2)

GO:0003407 neural retina development
Input ID	Gene Symbol	H Gene ID	Synonyms	Orphanet	OMIM
ENSDARG00000005374	tubgcp4	TUBGCP4	76P|GCP‐4|GCP4|Grip76|MCCRP3	[2518] Autosomal recessive chorioretinopathy‐microcephaly	OMIM:609610
ENSDARG00000098080	gnl2	GNL2	HUMAUANTIG|NGP1|Ngp‐1|Nog2|Nug2		OMIM 609365
ENSDARG00000052348	smarca5	SMARCA5	ISWI|SNF2H|WCRF135|hISWI|hSNF2H	[370334] Extraskeletal Ewing sarcoma	OMIM:603375
ENSDARG00000002235	mmp14a	MMP14	MMP‐14|MMP‐X1|MT‐MMP|MT‐MMP 1|MT1‐MMP|MT1MMP|MTMMP1|WNCHRS	[85196] Nodulosis‐arthropathy‐osteolysis syndrome;[3460] TORG‐WINCHESTER SYNDROME	OMIM:600754
ENSDARG00000018492	znf503	ZNF503	NOLZ‐1|NOLZ1|Nlz2		OMIM:613902
ENSDARG00000025187	six6a	SIX6	MCOPCT2|ODRMD|OPTX2|Six9	[264200] 14q22q23 microdeletion syndrome;[435930] Colobomatous optic disc‐macular atrophy‐chorioretinopathy syndrome;[2542] Isolated anophthalmia‐microphthalmia	OMIM:606 326
ENSDARG00000029688	hsp70.1	HSPA1L	HSP70‐1L|HSP70‐HOM|HSP70T|hum70t		OMIM:140559
ENSDARG00000069552	atoh7	ATOH7	Math5|NCRNA|PHPVAR|RNANC|bHLHa13	[289499] Congenital cataract microcornea with corneal opacity;[91495] Persistent hyperplastic primary vitreous	OMIM:609 875
ENSDARG00000071684	rx1	RAX	RX |MCOP3	[2542] Isolated microphthalmia‐anophthalmia‐coloboma	OMIM:601881

We next investigated the known relationships among the 137 Atoh7‐regulated genes via network reconstruction analysis (see methods). Interestingly, analysis conducted with STRING‐DB v 11.0  highlighted atoh7, rx1, and six6a as the core of a retinal “kernel” composed of early developmental eye‐specific transcription factors (Figure 4). Furthermore, this analysis highlighted three main gene subnetworks, which are suggestive of the over‐represented GO Biological processes, namely “retinal development,” “cell cycle/chromatin remodeling,” and “Wnt signaling pathway” (Figure 4).

FIGURE 4

Interaction network downstream of Atoh7. Known interactions between downstream targets of Atoh7 from the STRING database visualized with Cytoscape (genes without known interactions are not represented). Node colours represent the log2 fold‐change of gene expression in lakritz vs wild‐type eyes. Node borders are coloured by gene ontology annotation: “neural retina development” (green), “cell cycle process” (yellow), “Wnt signaling pathway” (pink), and “chromatin remodeling” (orange) [Color figure can be viewed at wileyonlinelibrary.com]

Among the bona fide genes integrating into the “cell cycle/chromatin remodeling” subnetwork is the F‐actin binding and cytokinesis regulator Anillin (ANLN , , ). Studies in the zebrafish retina have shown that anillin expression levels are required to favor cell cycle progression and restrict RGC genesis in Atoh7‐expressing RPCs. Concordantly, in the presence of anillin downregulation, many more RPCs turn on atoh7 and become RGCs, In agreement with these findings, our results show that anillin is an Atoh7‐downregulated gene (Figure 2). Studies have reported anillin expression positively associated with the expression of β‐catenin (ctnnb1). Interestingly, our results indicate that ctnnb1 is an Atoh7‐downregulated gene as well as a central gene of the “Wnt signaling pathway” emerging in the network reconstruction analysis (Figure 4). Given the reported role for Wnt/β‐catenin signaling in retinal development and regeneration, we sought to uncover possible interplay between Anillin and β‐catenin in the developing zebrafish retina. In support of this hypothesis, we find that anillin knockdown causes accumulation and displacement of the β‐catenin signal in the apical and apical‐lateral membrane of RPCs at 30 hpf (Figure 5A,B). Further studies linked anillin to angiogenesis and retinal neovascularization processes, leading to the obvious hypothesis that anillin regulation is relevant for retinal vascular development. If this were the case, we should find anillin expression being regulated beyond the RGC differentiation period. To begin to assess this possibility, we carried out quantitative real‐time RT‐PCR (qRT‐PCR) analysis of anillin starting from 25 hpf until after 72 hpf, when all retinal cells are fully differentiated in the central zebrafish retina. Interestingly, anillin expression levels remain significantly upregulated by the lakritz condition (Figure 5C) even at 72 hpf. Furthermore, the trend of anillin upregulation in subsequent developmental stages suggests oscillatory dynamics of anillin regulation during eye development (Figure 5C).

FIGURE 5

Dysregulation of Ctnnb1 localization by anillin knockdown and anillin expression dynamics. A, Ctnnb1 staining in control (ctrlMO, n = 3 embryos) versus anillin knockdown (anlnMO, n = 3 embryos) in morpholino injected embryos at 30hpf. Arrows show apical location, arrowheads show apical‐to‐basal location. B, Graph showing the normalized Ctnnb1 intensity signal along the basal‐to‐apical membrane of the apical‐most cells in control (n = 3 embryos) versus anlnMO (n = 3 embryos) injected embryos. The colored line shows the averaged intensity of all lines for the ctrlMO and anlnMO. Boxplot showing the number of peaks of Ctnnb1 signal intensity along the apical surface in ctrlMO (n = 3 embryos) versus anlnMO (n = 3 embryos) injected embryos. P < 10−4. Center lines show the medians; crosses show the means; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, data points are represented as circles. Student's t test. C, Anillin mRNA levels show dynamic variations during subsequent developmental stages. qRT‐PCR was performed on eyes from lakritz or wild‐type embryos at 25, 35, 48, 72 (left) and 96 hpf (right) to assess the trend of anillin and atoh7 expression. The relative gene expression (lakritz vs wild type) was calculated using the CT method for each stage. Histogram values are expressed as mean ± SEM. (P < .05) and the mRNA levels of both gapdh and ube2a were used as internal controls. The statistical analysis is described in the methods section. [Color figure can be viewed at wileyonlinelibrary.com]

Weighted Gene Co‐expression Network Analysis revealed a coexpression module with a cluster of genes tightly interconnected to Atoh7

To explore global changes of gene expression in the lakritz mutant, a weight gene co‐expression network analysis (WGCNA) was performed. Notwithstanding the small number of array samples, we identified 16 recurrent functional modules based on co‐expression pattern analysis on the full transcriptome dataset. To identify co‐expression modules significantly affected by the lakritz mutation, we tested their association with available covariates, including batch and mutation status. Out of the 16 modules found, two of them showed a high correlation with the mutation condition (wild type vs lakritz). Specifically, modules 13 (overall upregulated) and module 3 (overall downregulated) were found highly significant with a correlation of 0.97 (P = .002) and −0.99 (P = 1e‐5), respectively (Figure S2 and Table S4). Functional network analysis by STRING DB https://www.string‐db.org/) on M13 ‐ the smallest module, which we also found to contain atoh7 ‐ revealed an enrichment in “eye morphogenesis (blue color in Figure S3), “retina layer formation” (red color in Figure S3), and a cluster of genes previously found to be dysregulated in “light responsive, circadian rhythm processes” (PMID: 20830285 and PMID: 21390203) (light and dark green color in Figure S3).

Lastly, we analyzed in detail the topology of the Module 13. We identified four submodules with high within‐community connectivity, which show decreasing degree of connectivity from left to right (Figure 6). The submodule containing atoh7 (left) shows a densely interconnected cluster of genes with high topological overlap. These genes likely participate in common regulatory and signaling circuits including retina layer formation (eg, atoh7, rx1) and Wnt/β‐catenin signaling pathway (eg, fxd7a, tcf7l1b).

FIGURE 6

Detailed topology of Module 13. Each node represents a gene while a connection represents a co‐expression between two genes (only the first 200 edges in order of co‐expression weight were retained for visualization purposes). Submodules are shown with decreasing degree of connectivity from left to right. Highlighted edges represent the connection between retina layer formation and Wnt‐related genes [Color figure can be viewed at wileyonlinelibrary.com]

DISCUSSION

Our differential gene expression analysis of transcriptome data revealed 137 genes that are significantly differentially expressed between lakritz and wild‐type eyes from embryos at a developmental time‐point corresponding to the onset of RGC differentiation. We also applied multiple bioinformatics pipelines to perform a functional classification and network reconstruction of the differentially expressed genes. Notably, all methods here applied consistently highlighted “neural retina development” (GO:0003407) as the most biological pathway differentially affected by the lakritz mutation. Likewise, the interplay atoh7, rx1, and six6a—early developmentally regulated EFTs—consistently emerged as the “kernel” network of this cluster

The homeobox transcription factor Rx1 is well known for its evolutionarily conserved role in the generation and maintenance of multipotent RPCs during morphogenesis and differentiation of the vertebrate eye. , , , , , Mutant variants of RAX family genes have been linked to congenital developmental eye disorders, particularly microphthalmia, further confirming an early requirement for RPC proliferation and stemness. , , , Based on these findings, the data from our analysis suggest that balancing RPCs competence and RGC fate commitment requires rx1 downregulation by Atoh7. Notably, Atoh7‐mediated downregulation of rx1 is likely direct, since previous in silico analyses highlighted the presence of an Atoh7‐binding motive in the rx1 gene cis‐regulatory regions.

Conversely, six6 genes have been reported as RPCs competence factors, on the one hand suppressing stemness and proliferation via Wnt/β‐catenin signaling downregulation, on the other hand promoting expression of RGC differentiation genes. These findings are in agreement with the present study reporting six6a being an Atoh7‐upregulated gene. Remarkably, variant forms of SIX6 (SIX9/OPTX2) have been linked to congenital microphthalmia as well as to the development of glaucoma. , , , , , , , , Given that mutations in ATOH7 have been also associated with similar global eye disorders, , these observations strongly point at the importance to further our understanding of the interplay ATOH7, RX1, and SIX6 during eye development as well as to assess how disruption of this evolutionarily conserved genetic network might be linked to such eye disorders.

Besides atoh7, rx1, and six6a, significant differentially expressed genes annotated with “neural retina development” were tubgcp4, gnl2, smarca5, mmp14a, znf503, and hsp70.1. The Atoh7‐upregulated gamma‐tubulin complex protein 4 encoding gene tubgcp4 is of great interest. Variants of TUBGCP4 have been linked to autosomal‐recessive microcephaly with chorioretinopathy, which comprise a spectrum of eye developmental anomalies including microphthalmia, optic nerve hypoplasia, retinal folds, and absence of retinal vasculature. This raises the possibility that regulation of this gene might link Atoh7 to retinal‐vascular development as well as retinal neurogenesis. Genetic evidences for the GTPase and zinc finger transcriptional repressor encoding genes, gnl2 and znf503 (NOLZ1) as eye disorders‐related genes are still missing but studies support their functional requirement for retinal developmental processes, including proper cell cycle exit of RPCs during RGC differentiation. , Notably, tubgcp4, gnl2 and znf503 were implicated in the regulation of cytokinesis late in mitosis. , , This observation further indicates Atoh7 requirements for the regulation of cell cycle progression, at least in RPCs. The SWI/SNF chromatin remodeling factors have been reported as crucial regulators of the transition from multipotent to committed progenitor and differentiated cell states in multiple eye tissues, with potential implications for eye disorders. , , , The finding of smarca5 as an upregulated gene in the lakritz indicates the importance to address functional implications of this chromatin remodeling gene for atoh7‐related eye disorders. Furthermore, the “neural retina development” Gene Ontology category encompassed the reportedly stress response genes hsp70.1(HSPA1L) and mmp14 as significantly downregulated and upregulated, respectively, by the lakritz mutation. The crystallin related, heat shock Hsp70 family proteins are emerging as important regulators of RGC survival and regeneration as well as retinal vascular remodeling factors. , , , The intriguing finding that hsp70.1 is highly enriched among the Atoh7‐upregulated genes in the “neural retina development” gene cluster supports the idea that upregulation of these stress‐response proteins might be relevant also during RGC development. Lastly, the extracellular matrix remodeling factor Mmp14 is reportedly a crucial regulator of cell stemness and vascular remodeling. Studies have also implicated Mmp14a in retinal developmental processes such as RGC axon guidance and innervation of the optic tectum. , , Future studies will assess the functional implication of Atoh7‐mediated downregulation of mmp14a during vascular‐retinal development. Interestingly, our network analysis reveals that mmp14a is linked with known components of the Wnt signaling pathway, further underscoring the importance of the interplay of this pathway and Atoh7.

The “Wnt signaling pathway” is the second Atoh7‐dependent subnetwork emerging in our functional network analysis, which is centered around the ctnnb1 (β‐catenin) gene (Figure 4). Studies have shown that the Wnt/β‐catenin pathway promotes RPC proliferation and stemness while suppressing atoh7 activation and RGC differentiation. , , We here find that the main differentially expressed components of this pathway, namely ctnnb1, fxd7a, and tcf7l1b, are upregulated in the lakritz mutant, suggesting Atoh7 requirement in the Wnt/β‐catenin pathway downregulation (Figure 4). Concordantly, studies have reported that atoh7‐expressing RPCs contain low levels of expression in Wnt/β‐catenin pathway components when compared with non‐atoh7‐RPCs. This further suggestive of a negative feedback regulatory loop integrating Atoh7 and Wnt/β‐catenin signaling. Conversely, the planar cell polarity (PCP) signaling component celsr1a/flamingo, which has been reported as key regulator of neuronal cell differentiation, neurite outgrowth, and axon guidance, emerges as an Atoh7‐upregulated gene in our cohort. A number of studies reported dysregulation of Wnt signaling being associated with retinal diseases. , , , Likewise, Wnt, Fzd7/β‐catenin pathway has been reported as an important modulator of retinal vascular remodeling. Further research will clarify the genetic regulatory networks integrating Atoh7 and Wnt/β‐catenin signaling in controlling multiple eye tissue development in the vertebrate, and how their dysregulation might results in multiple ocular disorders. ,

Previous studies highlighted the importance of the interplay of Atoh7 with components of the Notch signaling pathway for RGC development and regeneration. , , In particular, downregulation of Notch signaling pathway has been proposed as a general mechanism whereby RGC genesis can be enhanced. , , , In agreement with these findings, Notch pathway components (hes6, mib1, hdac, adam17b, and Notch1a) appear affected by the lakritz mutation in our gene expression microarray data but their expression does not change significantly between wild‐type and lakritz condition in our analysis (Table S1). However, even though they fail short of the significance threshold of adjusted P value < .05, it is worth noting that their values suggest a trend in downregulation of Notch pathway‐related genes (Table S1). We suppose that such lack of significance in the differential expression might be linked to the developmental stage used for this analysis. Likewise, the extracted cohort of significantly regulated genes does not comprehend some of the reported RGC maturation‐associated factors, such as Cxcr4b, Elavl3, and Isl1. , , , Nonetheless, for cxcr4b and elavl3a, −1.43 FC (nominal P value of .0018) and −1.39 FC (nominal P value .0066) was observed, respectively, consistently with a positive regulation by Atoh7 (downregulated in the lakritz). We therefore suppose that expression of such RGC maturation‐related factors downstream of Atoh7 might be too low at the developmental time‐point selected for this study, to be detected within the chosen significance range. In support of the Atoh7‐dependent transcriptional regulation of RPC division and developmental progression, a third Atoh7‐dependent subnetwork emerged in our functional network analyses, which includes cell‐cycle and chromatin regulators (Figure 4). One bona fide gene of great interest in this subnetwork is the F‐actin binding and cytokinesis regulator Anillin (ANLN). , , Anillin has attracted increasing attention as a potential disease‐related gene (ORPH:93213, OMIM:616027). Evidence point at anillin expression levels being associated with cell proliferation and migration disorders in cancer and kidney diseases. , , , Additional roles for this actin binding protein have been reported in nerve cell development , and dysregulation of Anillin has been implicated in central nervous system myelin disorders. We here find anillin as an Atoh7‐downregulated gene. This is in line with earlier observations supporting the suggestion of a molecular feedback regulatory loop of an as yet unknown nature between anillin and atoh7 balancing RPCs developmental progression. While further investigations will address this question, the functional and network analyses from this study suggest that the Wnt/β‐catenin signaling might be involved. We indeed made preliminary observations indicating that anillin is required for the accumulation and/or localization of β‐catenin in the apical and apical‐lateral membrane of RPCs. Notably, in addition to promoting cell‐to‐cell adhesion , β‐catenin functions as nuclear transcriptional co‐activator of Wnt signaling responsive genes promoting proliferation and inhibiting differentiation. , , Thus, regulation of β‐catenin accumulation and localization to the E‐Cadherin/β‐catenin complex might be a mechanism whereby Anillin controls not only cell adhesion dynamics but also Wnt signaling pathway activity. It was also intriguing to note that anillin expression levels are affected by the lakritz mutation even after neurogenesis is completed in most of the central retina of the zebrafish eye. Given that Anillin has been reported as potential regulator of choroidal angiogenesis, , it is tempting to speculate that both RPCs and endothelial cell behaviors require anillin expression levels during retinal‐vascular developmental interactions. Future studies will assess the functional implications of the interplay of Atoh7, RGC genesis, Anillin and Wnt pathway components for the dysregulation of developmental vascular‐retinal disorders.

Finally, we have applied weighted gene co‐expression network analysis to explore gene co‐expression relationships and identify co‐expression modules potentially involved in atoh7 function. In addition to highlighting the already known interaction networks that were enriched within the 137 differentially regulated genes, we were able to identify two co‐expression modules significantly affected by the lakritz mutation (Table S4). One of them in particular (module 13) contains a cluster of highly interconnected genes including atoh7 itself, rx1 and members of the Wnt signaling pathway (eg, tcf7l1b, fzd7a and mmp14a). This analysis therefore further confirms tight functional interaction between neural retinal development and Wnt pathway genes from our gene differential expression data. Lastly, our knowledge‐based analysis of the M13 members by STRING database further extended these finding by revealing the existence of a network of interactions between atoh7‐ and rx1‐dependent “retina layer formation,” “eye morphogenesis” pathways, Wnt signaling pathway components (eg, tcf7l1b, fzd7a, and mmp14a), and two new gene networks previously found to be dysregulated in light responsive, circadian rhythm processes , (Figure S3). At present, very little is known on the functional importance of circadian clock genes, but increasing evidence indicates their implication in multiple eye tissues developmental processes as well as ocular disorders. , , , , This study further supports this evidence, by showing that Atoh7‐dependent regulatory networks integrates such circadian clock genes.

In sum, this Atoh7 targets analysis extends data from other studies focusing on transcription factors cascades enhancing RGC differentiation. We here provide new insights on Atoh7‐dependent developmental processes that might be regulated in global developmental eye disorders. First, they suggest that Atoh7 directly controls a two‐tiered regulatory network balancing early acquisition of progenitor cell competence (eg, through six6a, rx1) and repression of pro‐multipotency and proliferative processes (eg, through chromatin remodeling, cell cycle, and Wnt pathway regulation). Second, these data highlight for the first time many previously unreported cytoskeletal proteins, chromatin remodeling factors, stress‐response proteins, and even circadian clock genes as Atoh7‐regulated genes. Third, this analysis underscores both direct and potential functional genetic links of many of these factors to eye developmental disorders. This study thus contributes to laying the groundwork for the identification of key candidate molecules and their networks as potential targets for early eye disease detection and therapeutic applications.

AUTHORS’ CONTRIBUTIONS

L. Poggi designed research; G. Covello, A‐L Duchemin, FB Tremonti, L. Poggi and J. Ngai performed experiments; FJ Rossello, M. Filosi, F. Gajardo, E. Domenici, M. Eichenlaub, M. Ramialison performed bioinformatic analysis with inputs from JM Polo, D. Powell, M. L. Allende; L. Poggi and M. Ramialison wrote the manuscript with inputs from all authors.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Supplementary Material

Click here for additional data file.