Genes Involved in Drosophila melanogaster Ovarian Function Are Highly Conserved Throughout Evolution.

This work presents a systematic approach to study the conservation of genes between fruit flies and mammals. We have listed 971 Drosophila genes involved in female reproduction at the ovarian level and systematically looked for orthologs in the Ciona, zebrafish, coelacanth, lizard, chicken, and mouse. Depending on the species, the percentage of these Drosophila genes with at least one ortholog varies between 69% and 78%. In comparison, only 42% of all the Drosophila genes have an ortholog in the mouse genome (P < 0.0001), suggesting a dramatically higher evolutionary conservation of ovarian genes. The 177 Drosophila genes that have no ortholog in mice and other vertebrates correspond to genes that are involved in mechanisms of oogenesis that are specific to the fruit fly or the insects. Among 759 genes with at least one ortholog in the zebrafish, 73 have an expression enriched in the ovary in this species (RNA-seq data). Among 760 genes that have at least one ortholog in the mouse; 76 and 11 orthologs are reported to be preferentially and exclusively expressed in the mouse ovary, respectively (based on the UniGene expressed sequence tag database). Several of them are already known to play a key role in murine oogenesis and/or to be enriched in the mouse/zebrafish oocyte, whereas others have remained unreported. We have investigated, by RNA-seq and real-time quantitative PCR, the exclusive ovarian expression of 10 genes in fish and mammals. Overall, we have found several novel candidates potentially involved in mammalian oogenesis by an evolutionary approach and using the fruit fly as an animal model.

Introduction

Several examples of conservation of genes of the biochemical and molecular mechanisms involved in various biological functions have been previously documented, in particular between flies and mammals (McGary et al. 2010). For example, hox genes are involved in morphogenesis and cell differentiation in the fly and in the morphogenesis of the vertebrate embryo (Mark et al. 1997); the ovo gene is implicated in epidermal differentiation in the fruit fly and in hair formation in the mouse (Dai et al. 1998).

In Drosophila, contrary to mammals, oogenesis spans the entire life of a female, due to the presence of stem cells in the most anterior part of the ovary in a region called the germarium.

Although different at first sight, Drosophila and vertebrate oocyte formation and maturation might share more common features than previously thought (Matova and Cooley 2001). Firstly, meiosis is an important common step in the precursors of oocytes as well as packing the oocyte with RNA, proteins, and organelles that are crucial for the first steps of oogenesis. Moreover, within the animal kingdom, neo-oogenesis exists in several species, not only in insects but also in molluscs and teleost fishes (Hann 1927; Craig-Bennett 1931; Eggert 1931; Billard 1987; Wallace and Selman 1990; Grier et al. 2007). In Drosophila, many experiments have shown that the self-renewal of germline stem cells takes place in a niche, a physical structure made of signaling cells, that functions to maintain stem cell self-renewal and to prevent stem cell differentiation (Watt and Hogan 2000; Spradling et al. 2001). Renewal of the ovarian reserve during adult life is widely found in teleost fishes, and, in some salmonid species, there is an increase in the number of oogonia following spawning. In the rainbow trout (Oncorhynchus mykiss) postovulatory ovary, cell nests can be found in the germinal epithelium that contain oogonia, early diplotene oocytes and prefollicle cells (Grier et al. 2007). In medaka (Oryzias latipes), recently characterized ovarian cords within the germinal ovarian epithelium have been hypothesized to be reminiscent of the germarium from the Drosophila ovary (Nakamura et al. 2010). In mammalian ovaries, this type of niche has never been observed neither based on morphological studies nor molecular markers.

Recently it has been shown, by analyzing how primary oocytes are formed from germ cells in the mouse embryo that germ cell cysts exist in mouse ovaries similar to those found in Drosophila. Cells of these cysts are connected by ring canals and result from the synchronous divisions but incomplete cytokinesis of primordial germ cells. As in Drosophila, the future murine oocyte seems to be “nursed” by surrounding sister cells transferring some of their cytoplasm to this cell including Golgi, mitochondria, centrosomes and likely RNAs and proteins. Eventually, this unique cell will become a primary oocyte, whereas the other cells of the cyst die, like nurse cells in Drosophila (Pepling and Spradling 2001). In both Drosophila and mice, oocyte transcripts encode the so-called maternal genetic determinants that will support the earliest stages of embryogenesis (Laver et al. 2015). Mechanisms that control mRNA localization and regulate translation play central roles in specification of the anteroposterior and dorsoventral axes of the oocyte and future embryo and have been studied in detail in Drosophila (Lasko 2012).

Many genes have been shown to be involved in oogenesis in Drosophila, and some of them have an ortholog in the mouse and zebrafish that is also involved in oogenesis. It is the case with mos and Aurora, involved in meiosis, and the mRNA binding-proteins orb/cpeb1 (Zheng and Dean 2007; Extavour 2009; Schindler et al. 2012). We have previously identified genes encoding mRNA binding proteins that are highly expressed in the mouse oocyte, based on their phylogenetic relationship with fruit flies and nematode orthologs known to be involved in oogenesis in both species (Drouilhet et al. 2008). However, systematically identifying the orthologs of all fly genes in mammals and other vertebrates remains undone. In the present work, we have recovered a list of 971 fruit fly genes whose mutations affect ovary, egg, follicle cells, oocyte, and/or germarium according to the FlyBase database (flybase.org/) or reported in the literature to affect Drosophila oogenesis. We have then systematically looked for orthologs of these genes in the Ciona, zebrafish, coelacanth, lizard, chicken, and mouse to 1) investigate their phylogenetic evolution and the evolution of ovarian function from the fly to vertebrates, 2) identify genes that are specific of the fly and/or the insects as well as genes that are conserved between these species, and 3) identify putative orthologs of genes whose expression are enriched or exclusively expressed in the mouse and/or zebrafish ovary, using in silico and RNA-seq data, respectively. To further validate the evolutionary-conserved expression pattern of a subset of genes, we further investigated their expression using another model teleost species, the medaka, which is evolutionarily distant from the zebrafish.

Materials and Methods

In Silico Analyses

A list of genes of Drosophila melanogaster involved in female reproduction at the ovarian level was generated by compiling FlyBase lists of Drosophila mutants affecting ovarian phenotypes (http://flybase.org/). Five different affected tissues or cells were chosen to increment the list: Ovary, egg, follicle cells, oocyte, and germarium. An additional list of candidate genes described as being involved in oogenesis in the literature was also added (Spradling 1993; McKearin and Christerson 1994; Bastock and St Johnston 2008; Kumano 2012; Stein and Stevens 2014). The Drosophila gene list thus included 356 genes affecting the egg, 223 genes affecting the oocyte, 317 genes affecting follicle cells, 184 genes affecting the germarium, 688 genes affecting the ovary, and 337 genes involved in oogenesis; the latter have been previously reported in the literature. Once redundancy is taken into account, the fruit fly gene list was composed of 971 unique genes (supplementary table S1, Supplementary Material online).

Orthologs of D. melanogaster genes (genome version BDGP6) were extracted from the Ensembl database using the BioMart tool (http://www.ensembl.org/biomart/martview/) for chordates: Ciona intestinalis (genome version KH), Latimeria chalumnae (coelacanth, genome version LatCha1), Danio rerio (zebrafish, genome version GRCz10), Anolis carolinensis (anole lizard, genome version AnoCar2.0), Gallus gallus (chicken, genome version Gallus_gallus-5.0), and Mus musculus (mouse, genome version GRM38.p5) (supplementary table S1, Supplementary Material online). The orthology percentage of the 971 Drosophila genes list was also retrieved in the Ensembl metazoa database (http://metazoa.ensembl.org/biomart/martview/) for sponge (Amphimedon queenslandica, genome version Aqu1), cnidarians (Nematostella vectensis, genome version ASM20922v1), platyhelminthes (Schistosoma mansoni, genome version ASM23792v2), Ecdysozoa: Arthropods (Drosophila erecta, genome version dere_caf1; Anopheles gambiae, genome version AgamP4; Apis mellifera, genome version Amel_4.5), roundworms (Caenorhabditis elegans, genome version WBcell235), Lophotrocozoa: Molluscs (Octopus bimaculoides, genome version PRJNA270931), annelids (Capitella teleta, genome version Capitella teleta v1.0), and echinoderms (Strongylocentrotus purpuratus, genome version Spur_3.1).

For all mouse orthologs of Drosophila genes, expressed sequence tag (EST) expression data from all tissues were collected from the UniGene database (http://www.ncbi.nlm.nih.gov/unigene). The ratio of ovarian expression was calculated as the sum of ovary, oocyte, fertilized, and unfertilized ovum in transcripts per million divided by the sum of all adult female tissues. The ratio of testicular expression was calculated as the sum of testis, epididymis and prostate in transcripts per million divided by the sum of all adult male tissues (supplementary table S2, Supplementary Material online).

For zebrafish orthologs, RNA-seq data were extracted from the PhyloFish database (Pasquier et al. 2016). The ratio of ovarian expression was calculated as the ovary expression divided by the sum of all adult female tissues (supplementary table S3, Supplementary Material online). The ratio of testicular expression was calculated as the testis expression divided by the sum of all adult male tissues (supplementary table S3, Supplementary Material online). For medaka, orthologs were identified using the list of manually curated orthology relationships established between zebrafish and medaka, using the spotted gar as an outgroup of teleosts (Braasch et al. 2016). Characterization of the gene lists was performed using gene ontology (GO) and DAVID 6.8 (https://david.ncifcrf.gov/home.jsp) (da Huang et al. 2009a,b).

Biological Material

All experimental protocols were carried out in accordance with the European Directive 2010/63/EU on the protection of animals used for scientific purposes and the recommendations of the French Ministry of National Education, Higher Education and Research. The brain, heart, liver, spleen, skeletal muscle (quadriceps), kidney, bladder, uterus, and ovaries were collected from three adult female mice of the Swiss strain (Janvier, Le Genest St Isle, France). Testes were also collected from three adult male mice. Tissues were frozen in liquid nitrogen and stored at −80 °C until their use for RNA isolation.

For medaka and zebrafish, all sampled fishes originated from the INRA LPGP experimental facility. Fishes were reared and handled in strict accordance with French and European policies and guidelines of the INRA LPGP Institutional Animal Care and Use Committee (# 25M10), which approved this study. In medaka, all tissues were collected from eight adult females (0.431 ± 0.05 g) and eight adult males (0.425 ± 0.05 g). Zebrafish fish tissue originated from adult fish, as previously described (Braasch et al. 2016; Pasquier et al. 2017). All samples were immediately frozen in liquid nitrogen and subsequently stored at −80 °C until RNA extraction.

Gene Expression Analysis

Total RNA was extracted from murine tissues using TRI reagent, as recommended by the manufacturer’s instructions. Then, potential genomic DNA contamination was eliminated using RQ1 RNase-free DNase (Promega, Charbonnières-les-Bains, France), and the RNA quantity was assessed using a NanoDrop ND-1000 spectrophotometer (Nyxor Biotech, Paris, France). Reverse transcription (RT) was carried out on 1 µg of total RNA with the iScript Select CDNA synthesis kit (BioRad, Marnes La Coquette, France), according to the manufacturer’s protocol. Real-time PCR was performed in a total volume of 20 µl containing each primer at a final concentration of 125 nM (see sequences of specific primer pairs in table 1), 5 µl of RT reaction solution (diluted 1:25) and 10 µl of IQ SYBR Green Supermix (BioRad), according to the manufacturer’s instructions. Real-time PCR reactions were run on a MyiQ cycler (BioRad), with the following thermal conditions: Denaturation at 95 °C for 3 min, then 40 cycles of a three-step protocol (denaturation at 95 °C for 30 s, specific annealing at 60 °C for 30 s, and elongation at 72 °C for 30 s) followed by acquisition of a melting curve. Gene expression was normalized using the housekeeping gene Rpl19 with stable expression under our conditions. The normalized values of relative expression (R) were calculated according to the following equation: where Cq is the cycle threshold and E is PCR efficiency for each primer pair. The specificity of the amplified fragment was controlled by checking the amplicon size by electrophoresis. Normalized gene expression is expressed as mean ±SEM of three female/male mice.

Table 1

Mouse Oligonucleotide Sequences

Abbrev.	Accession Number	Gene	Forward Primer	Reverse Primer	bp	E %
Eif4e1b	NM_001033269	eukaryotic translation initiation factor 4E family member 1B	GGTTGGGAATTGTTTTGAGG	CCTGGTACCCAATGATGGTC	185	99.9
Foxr1	NM_001033469	forkhead box R1	CTTCCATCTGATTGCCTTGG	TAACGACCCAGGAGAATTGG	180	108
Kpna7	NM_001013774	Karyopherin alpha 7	TGGAGATGGTGACAGTTTGC	CACCAGAGACCATGTGGATG	173	103.9
Lin28a	NM_145833	lin-28 homolog A	TGGACACCAGAAAAGGGACT	GATGGGGAAAGGACAGACAA	211	102.7
Nanos2	NM_194064	Nanos homolog 2	CGACCAGGCTCATACACTCA	GAGCAGAGCAGAGGGGACTA	179	102.2
Pin1rt1	NM_001033768	Protein (peptidyl-prolyl cis/trans isomerase) NIMA-interacting 1, retrogene 1	CAAGGGGGTGACAACAGTCT	GCGTTGGTGATGTGATTGAA	193	95.8
Rnf17	NM_001033043	Ring finger protein 17	CCCTTAAGTGACTCGGGAAA	TCTGCAAGCCCCATTTCTAC	200	104.7
Rpl19	NM_009078	Ribosomal protein L19	CCTCCAGGCCAAGAAGGAAG	GGGCAACAGACAAAGGCTTG	151	87.8
Slc25a31	NM_178386	Solute carrier family 25 (mitochondrial carrier; adenine nucleotide translocator), member 31	GCCTTCCTAGAGACACCCTGT	CCCAAGCAACAGGAAAATGT	245	73.6
Spocd1	XM_017320505	SPOC domain containing 1	CAGGTGAAAATGGTCCTGCT	TCTTTGGGACCTTGTTCGAC	181	118.6
Tdrd1	NM_031387	Tudor domain containing 1	AAACAATCCAACCAACCAAAGT	GAAAGAGGAGGAAGCAGCAA	205	102.1
Tdrd9	NM_029056	Tudor domain containing 9	CTAGTGTGTGATGGGCCAAA	TCCACTCGTAGGGTTTCTCG	157	98.3

In medaka, total RNA from several tissues was extracted using the Trizol method (MRC, USA), according to the manufacturer’s instructions. The quantity and the quality of each RNA were assessed by measuring their absorbance at 260 and 280 nm using a Nanodrop 1000 spectrophotometer (ThermoScientific, Courtaboeuf, France) associated with ND-1000 V3 7.0 software. A DNAse treatment to remove residual DNA was carried out using the Ambion DNAse free kit, according to the manufacturer’s instructions. First-strand cDNA was synthesized with 540 ng of total RNA in each reaction system using the ThermoScientific Maxima first strand cDNA synthesis system, according to the manufacturer's protocol. At this step, luciferase control RNA (Promega), 50 pg per 500 ng of total RNA, was added to each sample to allow for data normalization. All cDNAs were diluted 1:10 with nuclease-free water and stored at −20 °C until used as templates in real-time quantitative PCR (RT-qPCR). Primers were designed so that they were overlapping an intron when it was possible (Primer3 software) using known sequences from the PhyloFish database. The primer sequences are shown in table 2. To determine the RT-qPCR efficiencies of each primer pair used, standard curves were generated using six serial dilutions (cDNA dilution from 1:5 to 1:160) of a pool from eight ovary samples. The ability of primers was validated when the amplification efficiency varied between 90% and 110%. Real-time quantitative PCR was carried out using a StepOnePlus Real-Time PCR System (Applied Biosystems, Villebon Sur Yvette, France) using Promega’s GoTaq qPCR Master Mix following the manufacturer's recommendations. Each PCR run included duplicates of each sample and also negative controls (reverse transcriptase-free and water RNA-free samples). The reactions were mixed in a volume of 10 µl containing 5 µl SYBR Premix, 4 µl cDNA (diluted 1:10) and 0.5 µl each of the 6 µM forward and reverse primers. After initial denaturation at 95 °C for 15 min, 40 cycles of amplification were carried out starting at 95 °C for 15 s, followed by 1 min at 60 °C, with a final extension at 60 °C for 1 min and 95 °C for 15 min. Data were subsequently normalized using exogenous luciferase transcript abundance in samples.

Table 2

Medaka Oligonucleotide Sequences

Abbrev.	Accession Number	Forward Primer	Reverse Primer
kpna7	XM_011494462.2	ACCAAAAGCCTTCTGTGCAG	GGGAAGGAGACCACATGTGA
bokb	XM_011486847.2	GGTGGTCCAAAGCCGAAATT	CACTCCAGCTCGTCTCCTAG
lsm14b-201	XM_011475365.2	AGTACCAGGAATTCAGGGCC	AACTCTGGATTTGGCGAGGA
mos	XM_020702078.1	GCACCTCCGTCATTCATTGG	CTCTTCACCTTCTTCAGCGC
pabpn1l-201	XM_011490104.2	ATCAAAGCCAGGGTTCAGGA	TAGAAAGGCCCAGGTTGAGG
cdt1-201	XM_004076871.3	ATTGTCTCAGACCACAGCCA	TGCTTCACGAGGGAGATCAG
orc1-201	XM_020702576.1	GGAAATCTGCGAGCACTCAG	CATTCAAGGCCTCCATCACG

Statistical Analysis

The proportion of genes exhibiting orthologs, enriched expression at the ovarian level, or enriched GO was compared using χ2 analyses.

Gene expression was compared between tissues using nonparametric analysis of variance (permutational ANOVA) (R package lmPerm) (Wheeler 2010) with the Tukey post hoc test (R package nparcomp; Konietschke et al. 2015), R version 3.3.1 (R Core Team 2015). A difference with P ≤ 0.05 was considered significant.

Results

Chordate Orthologs

The tables presenting the 971 fruit fly genes and their orthologs in chordate species, as well as EST/RNA-seq expression of corresponding murine/zebrafish ortholog genes are provided in supplementary tables S1–S3, Supplementary Material online. Among these 971 fruit fly genes, 81.8% (794 genes) possess a chordate ortholog in at least one of the following species: C. intestinalis (69.3%), L. chalumnae (coelacanth, 77.4%), D. rerio (zebrafish, 78.2%), A. carolinensis (anole lizard, 77.0%), G. gallus (chicken, 76.4%) and/or M. musculus (mouse, 78.3%). Considering only the mouse, 78.3% of the fruit fly genes (760 genes) have mouse orthologs, with a total of 1,556 mouse orthologs. Only 41.7% of the complete fruit fly gene repertoire (17,559) have mouse orthologs (7,318 genes) according to the Ensembl BioMart database. Therefore, our fruit fly gene list showed a significant increase in the proportion of genes possessing a mouse ortholog (P < 0.0001) or a vertebrate ortholog (fig. 1) in comparison to the fruit fly gene repertoire.

Fig. 1.

—Drosophila genes possessing a chordate orthologous gene. Number and percentage of orthologs of the 971 Drosophila melanogaster genes and of the complete Drosophila gene repertoire were obtained for each species considered using the Ensembl metazoa gene and Ensembl gene databases and the BioMart tool, as described in Materials and Methods section. Comparison of ortholog percentage of the 971 selected genes and the Drosophila gene repertoire were performed using χ2 tests. *Significant difference (P < 0.0001).

Given the high percentage of the studied genes that are conserved in the fruit fly and mouse, we investigated the conservation of these genes in other phyla. Among the 971 fruit fly genes, 66.0% possess an ortholog in sponge (A. queenslandica, 641 genes) and 72.0% possess an ortholog in cnidarians (N. vectensis, 699 genes). Only 34.6% and 38.6% of the complete fruit fly genome (17,559) have orthologs in A. queenslandica (6,067 genes) and N. vectensis (6,784 genes), respectively. Therefore, our fruit fly gene list showed a significant increase in the proportion of genes possessing an ortholog in sponge and cnidarians (P < 0.0001) compared with the fruit fly gene repertoire.

To determine which GO was especially enriched throughout evolution, a GO analysis was performed on the 1,556 mouse ortholog genes and the 1,849 zebrafish genes corresponding to the 971 fruit fly genes (figs. 2 and 3; supplementary tables S4 and S5, Supplementary Material online), and they were compared with the murine and zebrafish genome, respectively. Enriched GO terms corresponding to at least 20 of the 1,556 or 1,849 mouse or zebrafish genes, respectively, were investigated. Only GO terms with a Bonferroni corrected P value <0.01 were considered significantly enriched compare to the complete genome. In the mouse, concerning the biological process GO category, ontologies related to signal transduction (peptidyl-tyrosine phosphorylation, peptidyl-serine phosphorylation, regulation of protein phosphorylation, canonical Wnt signaling pathway, transmembrane receptor protein tyrosine kinase signaling pathway and cellular response to insulin stimulus) (P < 0.0001), regulation of translation (P < 0.0001) and to cell cycle (regulation of cell cycle and actin cytoskeleton organization) (P < 0.0001) were enriched more than 4-fold in the 1,556 mouse ortholog genes compared with the mouse genome (fig. 2; supplementary table S4, Supplementary Material online). Concerning molecular function GO category, ontologies related to kinase, endopeptidase, and phosphatase activities (P < 0.0001) as well as related to mRNA binding (P < 0.0001) were also enriched more than 4-fold in the 1,556 mouse gene list compared with the mouse genome (fig. 3; supplementary table S4, Supplementary Material online).

Fig. 2.

—Biological process GO enrichment. Fold enrichment of biological process GO in either zebrafish or mouse ortholog genes compared with the 971 Drosophila genes. GO analyses were performed using DAVID 6.8 (https://david.ncifcrf.gov/home.jsp). Only significantly enriched biological processes were represented (P < 0.05, Bonferroni corrected).

Fig. 3.

—Molecular function GO enrichment. Fold enrichment of molecular function GO in either zebrafish or mouse ortholog genes compared with the 971 Drosophila genes. GO analyses were performed using DAVID 6.8 (https://david.ncifcrf.gov/home.jsp). Only significantly enriched molecular functions were represented (P < 0.05, Bonferroni corrected).

In the zebrafish, similar ontologies were found significantly increased compared with the genome. Concerning biological process GO category, ontologies related to signal transduction (peptidyl-tyrosine autophosphorylation, integrin-mediated signaling pathway, canonical Wnt signaling pathway, transmembrane receptor protein tyrosine kinase signaling pathway, and peptidyl-serine phosphorylation) (P < 0.0001), cell migration (P < 0.0001), and actin cytoskeleton organization (P < 0.01) were enriched more than 3-fold in the 1,849 zebrafish ortholog genes compared with the zebrafish genome (fig. 2; supplementary table S5, Supplementary Material online). Concerning molecular function GO category, ontologies related to kinase and endopeptidase activities (P < 0.0001), motor activity (P < 0.0001), and microtubule binding (P < 0.0001) as well as related to mRNA binding (P < 0.0001) were also enriched more than 3-fold in zebrafish genes compared with the zebrafish genome (fig. 3; supplementary table S5, Supplementary Material online). A lot more GOs are enriched in the mouse gene list compared with the zebrafish gene list. Indeed, approximately all GOs that are enriched in the zebrafish gene list are also found enriched in the mouse gene list. Therefore, additional GOs are found enriched in the mouse gene list but not in the zebrafish gene list (i.e., regulation of translation, cellular response to insulin stimulus, positive regulation of neuron apoptotic process, negative regulation of ERK1 and ERK2 cascade, transforming growth factor β receptor signaling pathway, and histone deacetylase binding).

Drosophila Genes With No Orthologs in Chordate

Among the 971 fruit fly genes, 18.2% of the fruit fly genes have no chordate orthologs (177/971) (supplementary table S1, Supplementary Material online). We investigated the evolution of these genes in order to detail in the absence of these genes (fig. 4). Concerning the Ecdysozoa phylum, 94.4% of these genes have orthologs in the closest species of Drosophila (Drosophila erecta), which is a significantly higher conservation percentage compared with the complete D.melanogaster genome (72.8%). In other insects phyla (A. gambiae and A. mellifera, respectively), 57.9% and 52.8% of the 177 Drosophila genes without vertebrate orthologs, are present, which is not different from the average conservation of the complete D.melanogaster genome (51.8% and 47.4%, respectively). All the other invertebrate phyla analyzed showed significantly lower conservation percentages of the 177 genes list compared with the average conservation of the complete D.melanogaster genome: 12.4% have roundworm orthologs (C. elegans) versus 31.4%; concerning the Lophotrocozoa phylum, 23.6% of the genes have a mollusc ortholog (O. bimaculoides) versus 41.6% and concerning Echinoderms phylum, only 19.7% of the genes have a S. purpuratus ortholog versus 40.0%. This result could suggest that most of these genes are either lost after insect phyla or appeared throughout insect evolution.

Fig. 4.

—Drosophila genes possessing no chordate orthologous gene. Number and percentage of orthologs of the 178 Drosophila melanogaster genes exhibiting no ortholog genes in chordates were obtained for each species considered using the Ensembl metazoa gene and Ensembl gene databases and the BioMart tool, as described in Materials and Methods section. Comparison of ortholog percentages of the selected 178 genes and the Drosophila gene repertoire were performed using χ2 tests. *Significant difference (P < 0.0001).

In Silico Ovarian Enriched Expression in Mice

The UniGene database enabled us to class the mouse orthologs regarding their ovarian expression. Among the 1,556 mouse orthologs, 76 genes had an ovarian expression ratio above 0.5 (meaning that more than 50% of the total number of transcripts from all tissues of these genes are in the ovarian compartment) (supplementary tables S2 and S6, Supplementary Material online). Among these genes, 11 mouse genes are described as strictly expressed in the ovaries (Kpna7, Aurkc, Nanos2, Pabpn1l, T, Cts7, Eif4e1b, Tdrd1, Rbmxl2, Pin1rt1, and Spocd1), including 3 genes that are also strictly expressed in the testes of males (Kpna7, Tdrd1, and Rbmxl2).

Moreover, 52 genes had a testis expression above 0.5 (supplementary table S7, Supplementary Material online). Among these genes, 14 mouse genes are described as strictly expressed in the testis but not in the mouse ovary (Fscn3, Gm4064, Tdrd6, Hmgb4, Hus1b, 4921511C20Rik, Gm382, Enthd1, Meiob, Capza3, Gm732, Mycs, D1Pas1, and Rbm31y).

GO analysis was used to determine which GO was enriched when selecting genes on a preferential ovarian expression ratio above 0.5 (76 mouse genes) compared with the complete list of 1,556 mouse ortholog genes (supplementary table S8, Supplementary Material online). Enriched GO terms corresponding to at least 4 of the 76 mouse genes, respectively, were investigated. Only GO terms with a raw P value <0.05 were considered significantly enriched compare to the 1,556 mouse gene list. Concerning biological process GO category, the ontologies related to the meiotic cell cycle (P = 0.0026) and regulation of translation (P = 0.045) were significantly enriched (7.8- and 3.0-fold, respectively) compared with mouse genes (fig. 5). Concerning molecular function GO category, the mRNA binding ontology was significantly (P = 0.029) enriched (2.9-fold) in mouse genes, exhibiting a preferential ovarian expression compared with all 1,556 mouse genes.

Fig. 5.

—Biological process GO enrichment in ovarian genes. Fold enrichment of biological process GO in mouse ortholog genes exhibiting an ovarian expression compared with the complete list of mouse ortholog genes (1,556 genes). Fold enrichment of biological process GO in zebrafish ortholog genes exhibiting an ovarian expression compared with the complete list of zebrafish ortholog genes (1,849 genes). GO analyses were performed using DAVID 6.8 (https://david.ncifcrf.gov/home.jsp). Only significantly enriched biological processes were represented (P < 0.05, raw P values).

Ovarian Enriched Expression in Zebrafish (RNA-Seq)

The RNA-seq data obtained on zebrafish tissues allowed the classification of zebrafish orthologs based on their ovarian expression (supplementary table S3, Supplementary Material online). Among the 1,849 zebrafish orthologs, 73 genes had an ovarian expression ratio above 0.4 (meaning that more than 40% of the total number of transcripts from all tissues of these genes are in the ovarian compartment) (supplementary table S3, Supplementary Material online). As zebrafish data consisted of RNA-seq data and not EST data, a 0.4 ratio was chosen instead of 0.5. Indeed, contrary to UniGene data, expression in several tissues is observed for all the zebrafish genes studied, likely due to the efficiency of RNA-seq methodology. Interestingly, 20 genes were common between the mouse and zebrafish (supplementary table S6, Supplementary Material online). Indeed, 20 mouse genes and their 20 zebrafish orthologs had an ovarian preferential expression: Kpna7, Pabpn1l, Aurkb, Eif4e1b, Rnf17, Foxr1, Tdrd9, Pld6, Cpeb1b, Ccnb1, Kif23, Igf2bp1, Kif14, Ncapg, Lsm14b, Mybl2b, Mcm10, Mcph1, Dlgap5, and Cdc7.

GO analysis was used on genes showing a preferential ovarian expression ratio above 0.4 (73 zebrafish genes) compared with the complete list of 1,849 zebrafish ortholog genes (supplementary table S9, Supplementary Material online). Enriched GO terms corresponding to at least 4 of the 73 zebrafish genes, respectively, were investigated. Only GO terms with a raw P value <0.05 were considered significantly enriched compared with the 1,849 zebrafish gene list. Concerning biological process, the ontologies related to DNA methylation involved in gamete generation (P = 0.0005; 21-fold enriched), regulation of translation (P = 0.0028; 7.7-fold enriched), DNA replication (P = 0.011; 8.1-fold enriched), and microtubule-based movement (P = 0.011; 8.1-fold enriched) were significantly enriched compared with zebrafish genes (fig. 5). Concerning molecular function, the ontologies related to RNA binding (P = 0.004; 3.0-fold enriched), ATPase activity (P = 0.0038; 11.3-fold enriched) and microtubule motor activity (P = 0.018; 6.8-fold enriched) were significantly enriched compared with the 1,849 zebrafish genes.

Ontologies enriched in zebrafish ovarian genes and related to ATPase activity, microtubule-based movement, or microtubule motor activity were not enriched in mouse ovarian genes. On the contrary, ontologies enriched in mouse ovarian genes and related to the negative regulation of translation or gene silencing by RNA were not enriched in zebrafish ovarian genes.

In Vivo Confirmation of Ovarian Enriched Expression

Among the mouse orthologs exhibiting an in silico ovarian enriched expression, the expression of 11 genes was checked by qPCR: Eife1b, Foxr1, Kpna7, Lin28a, Nanos2, Pin1rt1, Rnf17, Slc25a31, Spocd1, Tdrd1, and Tdrd9 (fig. 6). Eight of these genes showed significantly enriched expression in the ovary compared with nongonadal tissues (Eife1b, Kpna7, Lin28a, Nanos2, Pin1rt1, Rnf17, Spocd1, and Tdrd1). Moreover, the expression in the ovary is higher compared with all other tissues for six of these genes: Eife1b, Kpna7, Lin28a, Nanos2, Pin1rt1, and Spocd1. Interestingly, a significantly enriched expression in the testis was found for six genes compared with nongonadal tissues: Kpna7, Rnf17, Slc25a31, Spocd1, Tdrd1, and Tdrd9.

Fig. 6.

—Gene expression in mouse tissues. Expression of 11 candidate genes was assessed in 10 mouse tissues (O, ovary; T, testis; M, muscle; U, uterus; Li, liver; S, spleen; B, brain; K, kidney; Bl, bladder and H, heart). After mRNA extraction and RT, qPCR was performed on tissues of three females or males. Expression of candidate genes is normalized by the expression of the Rpl19 housekeeping gene and presented as mean±SEM of three replicates. Different letters indicate a significant difference (P < 0.05). Therefore, if, at least, one letter is common between conditions, there is no significant difference between these two conditions (e.g., a is different from b but ab is not different from b. . .).

Among the genes exhibiting an ovarian enriched expression in zebrafish RNA-seq data, seven were selected for further analysis in medaka using qPCR: kpna7, pabpnl1, bokb, mos, lsm14b, cdt1, and orc1 (fig. 7). All genes had predominant expression in the ovary with the exception of bokb that was equally expressed in the eyes. Significant testicular expression was also observed for pabpnl1, cdt1 and orc1.

Fig. 7.

—Gene expression in medaka tissues. Expression of 7 candidate genes was assessed in 10 medaka tissues (O, ovary; T, testis; M, muscle; L, liver; B, brain; K, kidney; I, intestine; E, eye; G, gills and F, fins). After mRNA extraction and RT, qPCR was performed on tissues of eight females or males. Expression of candidate genes is normalized using exogenous luciferase transcript abundance in samples and presented as mean±SEM of eight replicates. Different letters indicate a significant difference (P < 0.05). Therefore, if, at least, one letter is common between conditions, there is no significant difference between these two conditions (e.g., a is different from b but ab is not different from b. . .).

Discussion

In the present study, we investigated the vertebrate orthologs of Drosophila genes involved in ovarian function. For the first time, we showed, from a list of 971 Drosophila genes, that around 80% of them were conserved throughout evolution. We also reported that 73 and 76 orthologs of these genes, in the zebrafish and mouse, respectively, were likely to exhibit preferential ovarian expression, 20 genes being common between both species.

Genes Shared by the Fly and the Vertebrates

Drosophila genes studied in this work showed high conservation throughout evolution. Indeed, an increased percentage of vertebrate orthologs of these Drosophila genes was found compared with the complete genome (80% versus 42%), in each species investigated. Moreover, orthology percentage was also increased in all other phyla investigated, including sponge, cnidarians, lophotrocozoa, ecdysozoa, and echinoderm, suggesting the critical importance of these genes. The other point to consider is that the orthology percentage we measured between vertebrate species and Drosophila is an underestimation. Indeed, due to the high stringency of the phylogenetic trees in Ensembl, we evidenced that already described mammal orthologs were not considered in the BioMart tool. For example, surprisingly, Vasa and Mos Drosophila genes were not predicted by BioMart to have an ortholog in vertebrates, although they do possess one (Gebauer and Richter 1997; Elis et al. 2008; Horie and Kotani 2016; Song et al. 2016). Interestingly, the conservation of “ovarian genes” is even higher for our ovarian genes than for “disease-focused genes” that averaged 65–70% from Drosophila to human (Hu et al. 2011; Wangler et al. 2017).

Among the GO enriched in the conserved genes, mRNA binding was reported as well as GO related to posttranslational regulation (phosphorylation) and signaling pathways. The processes of mRNA biogenesis, translation, alternative splicing, nuclear exportation, cytosolic localization, and finally degradation are modulated by RNA-binding proteins, which are of primary importance in the oocyte (Akagi et al. 2000). Several mechanisms contributing to this translational control, including conserved RNA-binding proteins that control translation by 3′ UTR binding and the assembly of mRNA–protein complexes called mRNPs, have been investigated not only in invertebrates but also in vertebrates (Drouilhet et al. 2008; Boateng et al. 2017; Pushpa et al. 2017; Rosario et al. 2017; Tsukamoto et al. 2017).

Conserved Genes Preferentially Expressed in the Ovary

In this study, among the fly genes that have orthologs, we have found that the expression of 73 and 76 of them are enriched in, or strictly specific to the ovary, in the zebrafish and mouse, respectively, including 20 genes that are preferentially expressed in the ovary in both vertebrate species. In the 73 and 76 zebrafish and mouse gene lists, respectively, the GO analysis suggested that mRNA binding function is increased compared with the full list of orthologs. Interestingly, among the 20 common genes between the zebrafish and mouse, 4 genes encode proteins involved in RNA binding and/or processing (pabpn1l, Lsm14, Eif4e1b, and Cpeb1). In particular, pabpn1l binds, stabilizes and localizes osk mRNA in Drosophila oocytes and might be involved in similar functions in vertebrates (Vazquez-Pianzola et al. 2011). The Lsm14 protein, known to associate with RNA to form the core domain of ribonucleoprotein particles, is involved in mRNA translation and distribution in Xenopus oocytes (Ladomery and Sommerville 2015). Eif4e1b was previously shown to be expressed in the mouse, Xenopus and zebrafish oocytes (Evsikov and Marin de Evsikova 2009) and to be involved in mRNA cap binding (Kubacka et al. 2015). Cpeb1 was also previously shown to regulate mRNA translation in the mouse (Sousa Martins et al. 2016) and Xenopus oocytes through the formation of a ribonucleic protein complex (Minshall et al. 2007). Four additional ovarian preferentially expressed genes, common between the zebrafish and mouse, can be related to RNA binding and processing (igf2bp1 and Rnf17, Tdrd9 and plD6). Indeed, the igf2bp1 gene encodes a protein that binds mRNA, but to our knowledge, its role in the oocyte is still not known. Rnf17, Tdrd9, and plD6 are all involved in piRNA processing. Indeed, the invalidation of Rnf17 (RING finger protein 17) leads to male sterility due to a blockage of spermiogenesis (Pan 2005), likely through the regulation of piRNA content and of PIWI proteins (Wasik et al. 2015). However, the same study showed that females are fertile. So, the role of this conserved gene in oogenesis remains to be determined. Moreover, Tdrd9 (tudor domain containing 9) encodes an ATP binding RNA helicase, and pld6 (or MITOPLD) is involved in the silencing of retrotransposons in male germ cells in the mouse (Shoji et al. 2009; Watanabe et al. 2011). A role of plD6 in piRNA processing was recently shown in the mouse oocyte, even if it seems to be dispensable in females, in contrast to males (Kabayama et al. 2017). In testes, molecular studies of germinal granules strongly implicate them in the transport, storage, localization, stability, and regulation of translation of mRNA (Houston and King 2000; Voronina et al. 2011).

Other functions are also represented in this common gene list between the zebrafish and mouse. Indeed, three genes encode serine/threonine kinases (AuroraC and Cdc7) and a cyclin (CCNB1). AuroraC and CCNB1 are known to play a role in the processes of mitosis and meiosis in oocytes (Uzbekova et al. 2008; Schindler et al. 2012), and Cdc7 is conserved from yeast to human and is essential for cell proliferation and embryonic development (Silva et al. 2006).

Six genes are involved in processes of cytokinesis, general chromosome relocation, DNA replication, recombination and repair and might play similar roles in the oocyte (Kif23, Kif14, Ncapg [or non-SMC condensin I complex, subunit G], Mcm10 [minichromosome maintenance 10], Mcph1 and Dlgap5 large [Drosophila] homolog-associated protein 5) (Geiman 2004; Liang et al. 2010; Lim et al. 2011; Liao et al. 2013; Liu et al. 2013; Samwer et al. 2013; Bonnet et al. 2015; Reubens et al. 2015). The gene encoding Kpna7, a karyopherin involved in the transport of molecules between cytoplasm and nucleus, has been shown to be necessary for the early steps of embryonic development in the pig (Wang et al. 2012). Lastly, two genes encode transcription factors: FOXR1, a forkhead DNA binding protein whose expression was described in human oocytes (Virant-Klun et al. 2013), but its role in oogenesis remains unknown; and Mybl2, involved in cell cycle progression and known to suppress autophagy in the pig ovary via the upregulation of the VDAC2 gene (Yuan et al. 2015).

Overall, these genes are highly conserved during evolution. Their specific ovarian expression seems to be well conserved as well and might play a universal role in oogenesis in both invertebrates and vertebrates. Of note, several of the other genes whose expression is highly enriched in the mouse and/or zebrafish ovaries, but under the threshold we fixed, have not been studied in the oocyte so far. Thus, future studies should also focus on their characterization in oocytes: Peptidyl-prolyl cis/trans isomerase Pin1rt1, Spocd1, T Brachyury transcription factor, Slc25a31, Naa11, or N(Alpha)-Acetyltransferase 11.

Limitations: False Negative Genes, No Ortholog Found

As mentioned before, one limitation of this study is that the stringency of the method used (Ensembl database and BioMart tool) prevented us, in some cases, from finding an ortholog even when we knew the vertebrate ortholog does exist. Therefore, it is possible that, concerning the supposed genes present in flies but not vertebrates, due to the evolutionary distance between these species, these genes have a true vertebrate ortholog, but that the high stringency of the phylogenetic trees recovered from the Ensembl database does not allow to identify them without any ambiguity (i.e., Vasa, mos, and brca2 are Drosophila genes known to have vertebrate orthologs). Moreover, for the same reason of evolutionary distance, the conserved synteny could not be used between Drosophila and vertebrates to complete or strengthen our phylogenetic trees.

Genes Present in the Fly and Not in Vertebrates

Nevertheless, in the list of 971 genes that are involved in ovarian function in the fly, 177 are supposed to have no ortholog in vertebrates. Some of these genes are known to be specific to oogenesis in flies, such as border cell formation and migration (six genes implicated in EGFR signaling, aos, grk. . . and four implicated in JAK/STAT signaling, apt, upd. . .), germ line stem cell functions (bam, ote. . .), oocyte fate determination (dap, ranshi, and stwl), ring canal formation with nurse cells (sosie, akap200. . .), dorsoventral polarity (spz, trk. . .), anteroposterior axis specification (bcd, osk, cup. . .), chorion and vitelline membrane formation (cp16, cp36, psd. . .) or dorsal appendage formation (Gfrl, brk. . .). One of the most interesting genes is oskar. Indeed in Drosophila, it seems that there is a hierarchical model for the organization of germ plasm, oskar being necessary and sufficient to initiate a cascade of interactions that assemble germ plasm at its site of localization (Mahowald 2001). Oskar initiates granule formation likely through direct interactions with Vasa, which in turn recruits Tudor (Breitwieser et al. 1996; Anne 2010), Valois, and PIWI factors (Kirino et al. 2009; Nishida et al. 2009). Germ granule mRNAs (nanos, germ cell less, and polar granule components) are transferred from the nurse cells to the oocyte (reviewed in Becalska and Gavis 2009). Surprisingly, some other genes are involved in more universal functions, such as apoptosis (hid, rpr, Strica. . .) or meiosis (mei-218, ord, cort. . .), suggesting that these genes are specific to these processes in Drosophila or insects.

Genes Involved in Oogenesis in the Fly and Expressed in the Testis in the Mouse

Surprisingly, some genes involved in oogenesis in the fly have a testis-specific expression in the mouse, suggesting a role in spermatogenesis, such as tra2/Rbmxl2 and tud/tdrd6. It is also the case in the ovo gene, which has been shown experimentally to be involved in spermatogenesis but not oogenesis in the mouse (Dai et al. 1998).

Conclusion

In conclusion, through this evolutionary in silico analysis using a D. melanogaster model, we have identified several vertebrate orthologs of genes involved in ovarian function and whose expression is highly enriched or specifically expressed in the ovary. Several of these genes are not known to play a role in oogenesis and represent new functional candidates to be tested by invalidation.

Funding

This work was supported by the IMAGE project (European H2020 Horizon framework programme, Grant agreement number 677353).

Click here for additional data file.