foxr1 is a novel maternal-effect gene in fish that is required for early embryonic success

Protein databases

Since our model is based on the zebrafish, all gene/protein nomenclatures will be based on those of fish. First, human (Homo sapiens) Fox protein R1 (Foxr1; ENSG00000176302) was used to BLAST for related protein sequences, and the following amino acid data were retrieved from the ENSEMBL database (http://www.ensembl.org/index.html): mouse, Mus musculus; rat, Rattus norvegicus; guinea pig, Cavia porcellus; pig, Sus scrofa; horse, Equus caballus; cow, Bos taurus; panda, Ailuropoda melanoleuca; opossum, Monodelphis domestica; Chinese softshell turtle, Pelodiscus sinensis; armadillo, Dasypus novemcinctus; frog, Xenopus tropicalis (Foxr1 a); fruit fly, Drosophila melanogaster; nematode, Caenorhabditis elegans; sea squirt, Ciona intestinalis; lamprey, Petromyzon marinus; coelacanth, Latimeria chalumnae; spotted gar, Lepisosteus oculatus (Foxr1 a); cod, Gadus morhua (Foxr1 a); fugu, Takifugu rubripes; medaka, Oryzias latipes (Foxr1 a); platyfish, Xiphophorus maculatus; stickleback, Gasterosteus aculeatus; tetraodon, Tetraodon nigroviridis; tilapia, Oreochromis niloticus; zebrafish, Danio rerio (Foxn1 and Foxn3); and cave fish, Astyanax mexicanus. Then, the protein sequence for zebrafish (Foxr1; NP_001096594.1) was found in the NCBI database (http://www.ncbi.nlm.nih.gov). Using this amino acid sequence as bait in BLAST, the bald eagle, Haliaeetus leucocephalu; penguin, Pygoscelis adeliae; crested ibis, Nipponia nippon; swan goose, Anser cygnoides domesticus; American alligator, Alligator mississippiensis; Chinese alligator, Alligator sinensis; python, Python bivittatus; central bearded dragon, Pogona vitticeps; frog, Xenopus laevis and Xenopus tropicalis (Foxr1 b) ; medaka, Oryzias latipes (Foxr1 b); northern pike, Esox lucius (Foxr1 a); rainbow trout, Oncorhynchus mykiss (Foxr1 a); coho salmon, Oncorhynchus kisutch; and Atlantic salmon, Salmo salar, protein sequences were extracted and investigated from the NCBI database. Further, using the protein sequence for zebrafish Foxr1, the following protein sequences were extracted from our previously established PhyloFish online database (http://phylofish.sigenae.org/index.html) (Pasquier et al., 2016) and analyzed along with the others: spotted gar, Lepisosteus oculatus (Foxr1 b); cod, Gadus morhua (Foxr1 b); bowfin, Amia calva; European eel, Anguilla anguilla; butterflyfish, Pantodon buchholzi; sweetfish, Plecoglossus altivelis; allis shad, Alosa alosa; arowana, Osteoglossum bicirrhosum; panga, Pangasius hypophthalmus; northern pike, Esox lucius (Foxr1 b); eastern mudminnow, Umbra pygmae; American whitefish, Coregonus clupeaformis; brook trout, Salvelinus fontinalis (Foxr1 a and b); rainbow trout, Oncorhynchus mykiss (Foxr1 b); European whitefish, Coregonus lavaretus; grayling, Thymallus thymallus; and European perch, Perca fluviatilis. These sequences are compiled in Data S1.

Phylogenetic analysis

The phylogenetic analysis was conducted using the Phylogeny.fr online program with default settings (Dereeper et al., 2008, 2010). Amino acid sequences of 73 Foxr1, Foxr2, Foxn1, and Foxn3 proteins from the above-mentioned species were aligned using the MUSCLE pipeline, alignment refinement was performed with Gblocks, and then the phylogenetic tree was generated using the Maximum Likelihood method (PhyML pipeline) with 100 bootstrap replicates.

Synteny analyses

Synteny maps of the conserved genomic regions of foxr1 and foxr2 were produced with spotted gar as the reference gene using PhyloView on the Genomicus v91.01 website (http://www.genomicus.biologie.ens.fr/genomicus-91.01/cgi-bin/search.pl).

Fish husbandry

Wildtype zebrafish (Danio rerio) of the AB strain were maintained at 25 °C in a central filtration recirculating system with a 12 h light/dark cycle in the INRA LPGP fish facility (Rennes, France). Individual couple pairing was performed by placing a male and a female overnight in a tank with a partition for separation, and in the morning, the divider was removed after which the female released her eggs to be fertilized by the male. All procedures of fish husbandry and sample collection were in accordance with the guidelines set by the French and European regulations on animal welfare. Protocols were approved by the Rennes ethical committee for animal research (CREEA) under approval no. R2012-JB-01.

Quantitative real-time PCR

Tissue samples from two wildtype males and three wildtype females, and 50–200 fertilized eggs at the one-cell stage from 32 wildtype couplings, which were all used as separate biological replicates, were harvested, and total RNA was extracted using Tri-Reagent (Molecular Research Center, Cincinnati, OH, USA) according to the manufacturer’s instructions. The tissue and egg samples were flash-frozen in liquid nitrogen immediately upon harvest and stored at −80 °C until use. Reverse transcription (RT) was performed using one μg of RNA from each sample with the Maxima First Strand cDNA Synthesis kit (Thermo Scientific, Waltham, MA, USA). Briefly, RNA was mixed with the kit reagents, and RT performed at 50 °C for 45 min followed by a 5-min termination step at 85 °C. Control reactions were run without reverse transcriptase and used as negative control in the quantitative real-time PCR (qPCR) study. qPCR experiments were performed with the Fast-SYBR GREEN fluorophore kit (Applied Biosystems, Foster City, CA, USA) as per the manufacturer’s instructions using 200 nM of each primer in order to keep PCR efficiency between 90% and 100%, and an Applied Biosystems StepOne Plus instrument. The StepOne software was used for expression analysis. RT products, including control reactions, were diluted 1/25, and four μL of each sample were used for each PCR. All qPCR experiments were performed in duplicate using technical replicates. The relative abundance of target cDNA was calculated from a standard curve of serially diluted pooled cDNA and normalized to 18S, β-actin, and EF1α transcripts. We considered a Ct variation of around 0.5 as acceptable. The R² values for foxr1, p21, p27, and rictor were 91.69%, 90.57%, 90.03%, and 90.36%, respectively. Melting curves were used and qPCR products were sequenced for confirmation. The primer sequences can be found in Data S2. The tissue expression of foxr1 was detected using the foxr1 genotyping forward and reverse primers while the mutant form of foxr1 in the CRISPR/cas9-mutated eggs was assessed with the foxr1 qPCR forward and reverse primers.

RNA-seq

Expression profiles in different holostean and teleostean species were obtained using the publicly available PhyloFish database http://phylofish.sigenae.org/index.html. Corresponding RNA-seq data were deposited into Sequence Read Archive of NCBI under accession references SRP044781–SRP044784, SRP045138, SRP045098–SRP045103, and SRP045140–SRP045146. The construction of sequencing libraries, data capture and processing, sequence assembly, mapping, and interpretation of read counts were all performed as previously described (Pasquier et al., 2016). The number of mapped reads was then normalized for the foxr1 gene across the 11 tissues using RPKM normalization.

In situ hybridization

Ovary samples were first fixed in 4% paraformaldehyde overnight, dehydrated by sequential methanol washes, paraffin-embedded, and sectioned to seven μm thickness before being subjected to the protocol. The sections were deparaffinized and incubated with 10 μg/mL of proteinase K for 8 min at room temperature, followed by blocking with the hybridization buffer (50% formamide, 50 μg/mL heparin, 100 μg/mL yeast tRNA, 1% Tween 20, and 5X saline-sodium citrate [SSC]). The probe was diluted to one ng/μL in the hybridization buffer and incubated overnight at 55 °C in a humidification chamber. The probes were synthesized by cloning a fragment of the foxr1 gene into the pCRII vector using the cloning foxr1 forward and reverse primers (Data S2) and Topo TA Cloning kit (Invitrogen, Carlsbad, CA, USA) as per the manufacturer’s protocol. The digoxigenin (DIG)-labeled sense and anti-sense probes were transcribed from Sp6 and T7 transcription sites, respectively, of the vector containing the cloned foxr1 fragment and purified using 2.5M LiCl solution. The purity and integrity of the probes were verified using the Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA) and the Agilent RNA 6000 Nano kit along with the Agilent 2100 bioanalyzer (Santa Clara, CA, USA). The slides were then subjected to two washes each with 50% formamide/2X SSC, 2X SSC, and 0.2X SSC at 55 °C followed by two washes with PBS at room temperature. The sections were subsequently blocked with blocking buffer (2% sheep serum, 3% bovine serum albumin, 0.2% Tween 20, and 0.2% Triton-X in PBS), and the anti-DIG antibody conjugated to alkaline phosphatase (Roche Diagnostics, Mannheim, Germany) was diluted by 1/500 and applied for 1.5 h at room temperature. The sections were washed with PBS and visualized with NBT/BCIP (nitro blue tetrazolium/5-bromo-4-chloro-3-indolyl phosphate).

CRISPR-cas9 genetic knockout

CRISPR/cas9 guide RNA (gRNA) were designed using the ZiFiT (Sander et al., 2010; Mali et al., 2013) online software and were made against two targets within the gene to generate a genomic deletion of approximately 240 base pairs (bp) that spans the last exon which allowed the formation of a non-functional protein. Nucleotide sequences containing the gRNA were ordered, annealed together, and cloned into the DR274 plasmid. In vitro transcription of the gRNA from the T7 initiation site was performed using the Maxiscript T7 kit (Applied Biosystems, Foster City, CA, USA) and of the cas9 mRNA using the mMESSAGE mMACHINE kit (Ambion/Thermo Scientific, Waltham, MA, USA) from the Sp6 site, and their purity and integrity were assessed using the Agilent RNA 6000 Nano Assay kit and 2100 Bioanalyzer. Zebrafish embryos at the one-cell stage were micro-injected with approximately 30–40 pg of each CRISPR/cas9 guide along with purified cas9 mRNA. The embryos were allowed to grow to adulthood, and genotyped using fin clip and PCR that detected the area around the deleted region. The full-length wildtype PCR band was 400 bp, and the mutant band with the CRISPR/cas9-generated deletion was approximately 160 bp using the foxr1 genotyping pair of primers. The PCR bands of the mutants were then sent for sequencing to verify the deletion. Once confirmed, the mutant females were mated with males harboring the vasa::gfp gene, where vasa was fused with gfp at the 3′ end, to produce F1 embryos, whose phenotypes were subsequently recorded. Images were captured with a Nikon AZ100 microscope and DS-Ri1 camera (Tokyo, Japan).

Genotyping by PCR

Fin clips were harvested from animals under anesthesia (0.1% phenoxyethanol) and lysed with 5% chelex containing 100 μg of proteinase K at 55 °C for 2 h and then 99 °C for 10 min. The extracted DNA was subjected to PCR using Jumpstart Taq polymerase (Sigma-Aldrich, St. Louis, MO, USA) and the foxr1 forward and reverse primers that are listed in Data S2.

Statistical analysis

Comparison of two groups was performed using the GraphPad Prism statistical software (La Jolla, CA, USA), and either the Student’s t-test or Mann–Whitney U-test was conducted depending on the normality of the groups based on the Anderson–Darling test. A p-value < 0.05 was considered as significant.

Phylogenetic analysis of Foxr1-related sequences

foxr1 and foxr2 were formerly known as foxn5 and foxn6, respectively (Katoh & Katoh, 2004a, 2004b). To date, there are six reported members of the foxr/foxn family (foxn1, foxn2, foxn3, foxn4, foxr1, and foxr2). To examine the evolution of foxr1, we used a Blast search approach using the zebrafish Foxr1 protein sequence as query in various public databases to retrieve 73 protein sequences from other species that could be related to this protein. All retrieved sequences are compiled in Data S1. Of note, both Foxr1 and Foxr2 protein sequences were retrieved. In order to verify that the retrieved protein sequences (Data S1) were homologous to zebrafish Foxr1, a phylogenetic analysis was performed. Based on the alignment of the retrieved vertebrate Foxr1-related sequences, and using Foxn1 and Foxn3 amino acid sequences as out-groups, a phylogenetic tree was generated (Fig. 1). As shown in Fig. 1, the common ancestor of the vertebrate foxr1 and foxr2 genes diverged from the ancestor of foxn1 and foxn3 genes, and these sequences were clearly observed as two separate clades belonging to actinopterygii (ray-finned fish) and sarcopterygii (lobe-finned fish and tetrapods). In addition, Foxr2 was found only in mammals with no homologs detected in actinopterygii as well as sauropsids and amphibians. Remarkably, despite the wide-ranging presence of Foxr1, no related sequences were observed in invertebrates and chondrichthyans (dogfish and sharks) or in certain species such as chicken (Gallus gallus). On the other hand, several species showed two Foxr1 sequences including the salmonids, rainbow trout (Oncorhynchus mykiss) and brook trout (Salvelinus fontinalis), as well as northern pike (Esox lucius), cod (Gadus morhua), medaka (Oryzias latipes), and spotted gar (Lepisosteus oculatus).

In line with the previous report that stated that foxr2 was absent in tilapia, stickleback, zebrafish, and medaka genomes, we retrieved Foxr2 protein sequences only inmammals using the zebrafish Foxr1 peptide sequence as query. Thus, using zebrafish Foxr1 sequence as the reference protein, we subsequently compared Foxr1 homology with the Foxr1 and Foxr2 sequences from mammals. As shown in Data S3, there was 29–37% identity and 41–53% similarity between all sequences, and there did not appear to be any difference in homology between zebrafish Foxr1 and mammalian Foxr1 and Foxr2 sequences. Further, there was 47–60% identity and 59–77% similarity between mammalian Foxr1 and Foxr2 sequences, indicating that these two proteins are highly similar and probably diverged recently during evolution.

Synteny analysis of foxr1 and foxr2 genes in vertebrates

In order to further understand the origin of the foxr1 and foxr2 genes in vertebrates, we performed a synteny analysis of their neighboring genes in representative vertebrate genomes using the basal actinopterygian, spotted gar, as the reference genome and the Genomicus online database (Fig. 2). We found that between the spotted gar and mammals, there was conserved synteny of the foxr1, upk2, ccdc84, rps25, trappc4, slc37a4, and ccdc153 loci in their genomes. In the frog (Xenopus tropicalis) genome, the foxr1, ccdc153, cbl, mcam, and c1qtnf5 loci were conserved, while in Coelacanth, foxr1, ccdc84, rps25, trappc4, slc37a4, cbl, ccdc153, mcam, c1qtnf5, as well as rnf26 loci were found in the same genomic region as those of the spotted gar. However, amongst the actinopterygians, there was lower conservation of synteny; in zebrafish and cave fish, the foxr1, ccdc84, and mcam loci were conserved while in the other ray-finned fish species, only the foxr1 loci was found. We further analyzed the foxr2 sequences that were found only in mammals, and we demonstrate here that they were all present on the X chromosome with no apparent conserved synteny of neighboring genes to those found in the spotted gar. Our overall analyses suggest that all the foxr-related sequences that were found were homologs, and the foxr gene in fish species probably derived from the ancestral foxr1 gene. Although there was the same degree of protein homology between zebrafish Foxr1 and mammalian Foxr1 and Foxr2 sequences, the phylogenetic tree and synteny analyses showed a clear distinction between them, and the foxr2 gene probably derived from a later single gene duplication or transposon event as previously suggested (Katoh & Katoh, 2004a).

Expression profiles of foxr1

We next focused our efforts on foxr1 since it has previously been shown in eel, tilapia, and medaka to be gonad specific and thus may have specific functions in reproduction. In order to investigate the potential functions of foxr1, we explored its tissue distribution using two different approaches, qPCR and RNA-seq, the latter of which was obtained from the PhyloFish online database (Pasquier et al., 2016). In zebrafish, we observed from both sets of data that foxr1 mRNA was predominantly expressed in the ovary and unfertilized egg (Figs. 3A and 3B). A high expression of foxr1 transcript was also detected in bone by RNA-seq, but not by qPCR; this could be due to methodological differences as RNA-seq is more sensitive than qPCR since it is designed to detect the sequences all along the transcript while the latter detects just one area of the transcript. Thus, it is possible that a different splice variant of foxr1 exists in bone. By in situ hybridization (ISH), we also demonstrated that foxr1 transcripts were highly expressed in the ovary in practically all stages of oogenesis (Figs. 3C–3E; negative controls, Figs. 3F–3H).

Functional analysis of foxr1 in zebrafish

To understand the role of foxr1 during oogenesis and early development, we performed functional analysis by genetic knockout using the CRISPR/cas9 system. One-cell staged embryos were injected with the CRISPR/cas9 guides that targeted foxr1 and allowed to grow to adulthood. Mosaic founder mutant females (F0) were identified by fin clip genotyping and subsequently mated with vasa::gfp males, and embryonic development of the F1 fertilized eggs was recorded. Since the mutagenesis efficiency of the CRISPR/cas9 system was very high, as previously described (Auer et al., 2014; Gagnon et al., 2014), the foxr1 gene was sufficiently knocked-out even in the mutant mosaic F0 females. This was evidenced by the substantially lower transcript level of foxr1 in the F1 embryos as compared to those from control pairings (0.78 ± 0.22 and 1.65 ± 0.04, respectively; Fig. 4A). Thus, the phenotypes of foxr1 (n = 5) mutants could be observed even in the F0 generation. Since neither males nor the females could transmit the mutated genes to future generations (i.e. all the surviving embryos were WT), all of our observations were obtained from the F1 generation, which were the fertilized eggs derived from the mosaic female foxr1 mutants.

We observed that most of the embryos from the foxr1 mutant females had a very low developmental success at 24 hpf (23.2 ± 8.4% vs. 85.2 ± 9.2% in controls; p < 0.008) (Fig. 4B). The frequency of the mutation in the mutant females is demonstrated in Fig. 4C, and it was observed that three of the mutants produced abundant non-developing eggs that remained non-cellularized, reflecting their failure to undergo cell division (Figs. 5E–5H). The eggs derived from these three foxr1 mutant females (Fig. 4: foxr1-1, foxr1-2, and fox1-3) did not undergo any cell division at two hpf and continued to display a complete lack of development up to eight hpf. By 24 hpf, these non-developing eggs that failed to divide were all dead. In addition, two of the mutants (Fig. 4: foxr1-4 and fox1-5) produced developmentally incompetent eggs with two phenotypes; those with a non-cellularized morphology (Figs. 5E–5H), and another population that developed albeit with an abnormal morphology (Figs. 5I–5L). These fertilized and developing embryos were structurally abnormal, with unsmooth and irregularly-shaped yolk as well as asymmetrical cell division that culminated into a blastodisc with a group of cells on top of an enlarged syncytium (Fig. 5K, arrow). These eggs underwent developmental arrest at around four hpf or the MBT and appeared to regress with further expansion of the syncytium (Figs. 5J–5K) until death by 24 hpf.

The observed phenotype of the foxr1 mutant-derived uncellularized eggs was very similar to previously described unfertilized eggs (Dekens et al., 2003). Thus, the foxR1 mutant females were mated with vasa::gfp males, and the genotype of their progeny was assessed for the presence of the gfp gene, which would only be transmitted from the father since the mutant females did not carry this gene. We found that these uncellularized eggs from the foxr1 mutant females did indeed carry the gfp gene (Fig. 5M), which indicated that some or all of them were fertilized, but were arrested from the earliest stage of development and did not undergo any cell division. These novel findings showed for the first time that foxr1 is essential for the developmental competence of zebrafish eggs, and is therefore a crucial maternal-effect gene.

In order to delve into the possible mechanisms that may be involved in the reduced reproductive success of the foxr1 mutants, we investigated the expression levels of p21, p27, and rictor, which were previously reported to be repressed by the Foxr1 transcription factor in mice (Santo et al., 2012). We found that there was substantially increased expression of p21 (4.82 ± 1.09 vs. 0.25 ± 0.04 in controls; p < 0.002) while that of rictor was significantly decreased (0.87 ± 0.13 vs. 1.60 ± 0.25 in controls; p < 0.01) in the foxr1 mutant-derived eggs as compared to eggs produced by wildtype females (Figs. 6A–6C). The expression of p27 was unchanged between the two groups. These results were in line with the growth arrested phenotype that was observed in the uncellularized and developmentally challenged eggs from the foxr1 mutant females.