Bo Zhu1
Wei Ge- 1Centre of Reproduction, Development and Aging, Faculty of Health Sciences, University of Macau, Taipa, China
- 2Genomics and Bioinformatics Core, Faculty of Health Sciences, University of Macau, Taipa, China
In teleosts, the onset of puberty in females is marked by the appearance of the first wave of pre-vitellogenic (PV) follicles from the pool of primary growth (PG) follicles (follicle activation) in the ovary during sexual maturation. To understand the mechanisms underlying follicle activation and therefore puberty onset, we undertook this transcriptomic study to investigate gene expression profiles in the event. Our analysis revealed a total of 2,027 up-regulated and 859 down-regulated genes during the PG-PV transition. Gene Ontology (GO) analysis showed that in addition to basic cellular functions such as gene transcription, cell differentiation, and cell migration, other biological processes such as steroidogenesis, cell signaling and angiogenesis were also enriched in up-regulated genes; by comparison, some processes were down-regulated including piRNA metabolism, gene silencing and proteolysis. Further Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis identified a variety of signaling pathways that might play pivotal roles in PG-PV transition, including MAPK, TGF-β, Hedgehog, FoxO, VEGF, Jak-STAT, and phosphatidylinositol signaling pathways. Other pathways of particular interest included endocytosis and glycosaminoglycan biosynthesis. We also analyzed expression changes of genes expressed in different compartments viz. oocytes and follicle cells. Interestingly, most oocyte-specific genes remained unchanged in expression during follicle activation whereas a great number of genes specifically expressed in the follicle cells showed significant changes in expression. Overall, this study reported a comprehensive analysis for genes, biological processes and pathways involved in follicle activation, which also marks female puberty onset in the zebrafish when occurring for the first time in sexual maturation.
Introduction
In each mammalian ovarian cycle, a limited number of primordial follicles are recruited to undergo a gonadotropin-independent growth phase (primary and secondary follicles) leading to antral formation (initial recruitment). This is followed by gonadotropin-dependent selection toward formation of pre-ovulatory follicles (cyclic recruitment) (1). As a critical stage of ovarian folliculogenesis, follicle recruitment or activation has been extensively studied in different mammalian models. Despite this, our understanding of the event and its controlling mechanisms is still limited (2), partly due to the limited supply of early follicles for study.
Similar to that in mammals, the process of folliculogenesis in fish has also been divided into gonadotropin-independent stage (primary growth, PG) and gonadotropin-dependent stage (secondary growth, SG), which are equivalent to preantral and antral stages in mammals respectively (3), although this view has recently been revised by emerging genetic data (4). During fish life cycle, the first transition from PG to SG in the ovary is considered an early visible marker for puberty onset or activation of the hypothalamic-pituitary-gonad (HPG) axis in females (5–7). Although the PG-SG transition has been described in many fish species, the mechanisms that control this important developmental event remain poorly understood. Studies in different fish species have provided lines of evidence for roles of pituitary hormones especially gonadotropins (FSH and LH), gonadal steroids and growth factors in regulating the event. A classical experiment of hypophysectomy in the goldfish showed an arrest of folliculogenesis at the transition from PG to previtellogenic (PV, the beginning of SG with cortical alveoli) stage without pituitary, strongly suggesting important roles for pituitary hormones in controlling follicle activation (8). This idea has been supported by studies in other species. For example, PG follicle growth and PG-PV transition in salmons were associated with increased plasma FSH level and gonadal FSH receptor expression (9, 10). Similar results have also been reported in the zebrafish (11, 12). Our recent study in the zebrafish using gene knockout approach provided further genetic evidence for the importance of gonadotropin signaling, especially FSH-FSHR pathway in controlling follicle activation (4, 13). In addition to pituitary gonadotropins, gonadal steroids especially androgens have also been implicated in regulating early follicle growth and possibly activation as well. In the shortfinned eel (Anguilla australis), treatment with fish-specific androgen 11-ketotestosterone (11-KT) significantly promoted the growth of PG follicles in vitro (14) and PG-PV transition in vivo (15). Similar results have also been reported in the Atlantic cod (Gadus morhua L.) (16, 17). Recent in vitro studies in the coho salmon (Oncorhynchus kisutch) confirmed the stimulatory effect of 11-KT on PG follicles and their transition to PV follicles. Interestingly, these studies provided evidence for sequential roles of androgens and estrogens in early follicle development. Androgens [11-KT and testosterone (T)] but not estrogens (E2) promoted the growth of PG follicles and their transition to PV stage; however, E2 exhibited strong stimulatory effect on the development of PV follicles after the transition (18, 19). Evidence has also accumulated for roles of growth factors in controlling early folliculogenesis. Both IGF-I and II increased their plasma levels during PG-PV transition together with FSH and E2 (10). In agreement with this, IGF-I, but not insulin and leptin, increased the size of eel PG follicles in vitro (14).
As a daily spawner with continual development of ovarian follicles, the zebrafish has become a popular model for studying follicle development and its controlling mechanism (20). In a recent study, we demonstrated that the female zebrafish entered puberty when the first cohort of PV follicles appeared in the ovary. This happened when the body weight reached 100 mg or standard body length above 1.8 cm, suggesting a role for the somatotrophic or growth axis in activating the reproductive axis (7), similar to the situation in salmons (10) and mammals (21). By using candidate gene approach, a variety of genes especially growth factors and their receptors have been shown to display differential expression patterns during the PG-PV transition in zebrafish, including activin-inhibin-follistatin family (22, 23), transforming growth factor β (TGF-β) family (24), bone morphogenetic protein (BMP) family (25–28), epidermal growth factor (EGF) family (29–31), insulin-like growth factor (IGF) family (32–34), Kit family (35), and anti-Müllerian hormone (AMH) (36). The candidate gene approach has recently been boosted by the emergence of powerful genome editing technologies such as TALEN and CRISPR/Cas9. An increasing number of genes have been analyzed for their functions in early folliculogenesis in fish models. For example, disruption of oocyte-specific transcription factor figla in zebrafish (37) and ovarian aromatase (cyp19a1a) in zebrafish (28, 38, 39), medaka (40) and tilapia (41) blocked early folliculogenesis and induced female-to-male sex reversal. In addition, inactivation of foxl2a/b in zebrafish (42) and tilapia (41), sf1 in tilapia (43), foxl3 in medaka (44), and bmp15 in zebrafish (28) all resulted in female-to-male sex reversal. Interestingly, the loss of fshr in female zebrafish caused complete failure of follicle activation or PG-PV transition followed by sex reversal to males (4, 45).
While the candidate target gene approach is still the major form of study for gene functions in folliculogenesis, there has been an increasing number of studies using high-throughput approaches such as transcriptomics using either microarray or next generation sequencing (NGS)-based RNA sequencing (RNA-seq) to investigate differentially expressed genes (DEGs) and regulatory networks. In humans, thousands of genes were shown by microarray analysis to be differentially expressed between mature oocytes and the somatic tissues (46). Similar approach also revealed hundreds of DEGs between the cumulus granulosa cells and the granulosa cells in the follicular fluid (47). Similar transcriptome studies have also been performed in a variety of mammalian species on temporal expression profiles at different stages of folliculogenesis, including bovine follicles of different sizes (48), primary and secondary follicles in sheep (49), secondary and early antral (tertiary) follicles in the goat (50), three different stages of antral follicles (small, medium and large) in the pig (51), and unassembled, primordial and primary follicles in the rat (2). In contrast to mammalian models, similar studies are limited in non-mammalian species. In the pigeon, pre-ovulatory and post-ovulatory ovaries were analyzed for transcriptomes, and more than 400 DEGs were identified (52). In teleosts, a few studies on ovarian transcriptomes have been reported using subtractive hybridization, serial analysis of gene expression (SAGE), microarray, or RNA-seq. In the flatfish, microarray analysis on three stages of follicle development, i.e., oocyte growth (vitellogenesis), maturation and follicle atresia, identified 118 DEGs (53). In the rainbow trout, three stages of precocious ovaries were analyzed for DEGs as compared to normal ovary by subtractive hybridization and microarray (54). A recent RNA-seq study in the common carp examined gene expression profiles at three stages of ovarian development, namely undifferentiated gonads, immature juvenile ovary and mature adult ovary, with a large number of DEGs assigned to biological processes related to reproduction (55). In the zebrafish, an early study used SAGE to analyze the transcriptome of full-grown follicles and compared it with similar transcriptomes of other vertebrate species (56). Despite these studies, there is a lack of transcriptome analysis on early folliculogenesis, which is critical for understanding the underlying mechanisms controlling follicle activation and therefore puberty onset. Using subtractive hybridization, a study in the coho salmon compared PG (perinucleolar stage) and PV follicles (early stage of SG with cortical alveoli) with the aim to identify DEGs between the two stages in early folliculogenesis. Some differentially expressed genes were identified, including zona pellucida glycoprotein (zp) genes, vitellogenin receptors (vldlr), anti-Müllerian hormone (amh), gonadal soma-derived growth factor (gsdf ), and follicle-stimulating hormone receptor (fshr) (57). Also in the coho salmon, a recent study demonstrated significant changes in ovarian transcriptome in response to 11-KT, which stimulates PG-PV transition. The DEGs identified included those involved in gonadotropin, steroid and growth factor signaling (19).
Using zebrafish as the model, we undertook this transcriptomic study to analyze gene expression profiles during follicle activation/recruitment or PG-PV transition by NGS-based RNA-seq. This high-throughput approach of gene expression analysis promises to provide us a comprehensive and global view on genes and regulatory networks or pathways involved in follicle recruitment or activation. The data obtained from the present study together with the availability of genome editing technologies in the zebrafish will make this species an excellent platform for studying gene functions in ovarian folliculogenesis.
Materials and Methods
Animals
The wild type zebrafish (Danio rerio) of AB strain was used in this study. The fish were maintained in the ZebTEC multilinking rack system (Tecniplast; Buguggiate, Italy) under an artificial photoperiod of 14 h light:10 h dark. The temperature, pH and conductivity of the system were 28 ± 1C°, 7.5 and 400 μS/cm, respectively. The fish were fed with paramecia, brine shrimp, and Otohime fish diet (Marubeni Nisshin Feed, Tokyo, Japan) twice a day during larval, juvenile and adult stages, respectively.
All experiments were performed under a license from the Government of the Macau Special Administrative Region (SAR), and approved by the Animal Experimentation Ethics Committee of the University of Macau.
Histological Examination
For histological examination, the dissected ovarian samples were immediately fixed in Bouin's solution overnight at room temperature. The fixed samples were dehydrated and embedded in paraffin according to our previous report (6), and then serially cut into 7-μm sections on a microtome (Leica, Wetzlar, Germany). Slides were stained with hematoxylin and eosin (H&E), followed by observation and photographing on Nikon ECLIPSE Ni-U microscope (Nikon, Tokyo, Japan).
Follicle Isolation and RNA Preparation
Zebrafish were anesthetized by ice shock and decapitated before dissection. Before follicle isolation, ovaries were carefully removed from three female zebrafish and placed in a 100-mm Petri dish containing 60% Leibovitz L-15 medium (Gibco Invitrogen, Carlsbad, CA). Fat tissue and ligaments surrounding the ovaries were stripped off using a BD Microlance 26G needle (BD, San Diego, CA). To separate follicles, a transfer pipette (JET Biofil, Guangzhou, China) was used to pipet ovarian fragments up and down a few times, followed by further pipetting with a BD 23G needle for a few times. After separation of follicles, the PG (<100 μm) and PV follicles (~180–250 μm) were isolated by filtering through sieves with different pore sizes, including 100 μm (SPL Lifesciences, Waunakee, WI), 180 and 250 μm (Jiu Feng, Heng Shui, China). The PG follicles that passed through the 100 μm sieve were collected by centrifugation at 3,000 rpm for 2 min, and the PV follicles retained between 180 and 250 μm were collected also by centrifugation. The PG and PV follicles were washed twice with 1 X PBS. In total, nine zebrafish were used and divided into three groups (3 fish per group) for follicle isolation and RNA preparation (n = 3). The entire process of ovary dissection, follicle isolation and washing was finished within 1–1.5 h to minimize changes of gene expression in vitro. Total RNA from PG and PV follicles of each group was extracted using Tri-Reagent (Molecular Research Center, Cincinnati, OH) according to the protocol of the manufacturer and our previous report (58). The RNA was then treated with DNase for 10 min at 37°C to remove genomic DNA [10 μg RNA in 100 μl reaction buffer with 2U DNase I from NEB (Ipswich, MA)] followed by phenol-chloroform extraction and ethanol precipitation.
RNA Sequencing and Differential Expression Analysis
For RNA library construction, the integrity of the RNA samples was first analyzed on the Bioanalyser 2100 (Agilent, Stockport, UK). RNA libraries were prepared using the NEBNext Ultra Directional RNA Library Prep Kit (NEB) and sequenced on the HiSeq 2500 Sequencing System (Illumina, San Diego, CA) with 100-bp paired end reads. FastQC (v0.11.5) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for quality checking of the raw fastq data. These raw reads were aligned against zebrafish reference genome (assembly GRCz10) using TopHat2 (v2.1.1) (https://ccb.jhu.edu/software/tophat/manual.shtml) (59) with default parameters. StringTie (v1.3.3b) (https://ccb.jhu.edu/software/stringtie/index.shtml) (60, 61) and prepDE.py Python script provided with StringTie tool were used to assemble the alignments into transcripts and extract the raw read counts for reference genomic features respectively. The read count matrix was processed by DeSeq2 package (v1.16.1) (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) (62) for differential gene expression analysis. Genes with p-values <= 0.05 and log2 (fold change) >= 1 were considered as differentially expressed genes (DEGs) with statistical significance.
The transcriptome data analysis was performed at the Genomics and Bioinformatics Core (Faculty of Health Sciences, University of Macau), and all the raw data were submitted to NCBI GEO database with accession number GSE99308.
Functional Enrichment Analysis
Significant DEGs were subject to GO and KEGG pathway enrichment analysis using online resource DAVID (v6.8) (63). GO terms and KEGG pathways with p-values < 0.05 were considered significantly enriched. DEGs with significant Log2 (fold change) were mapped onto the enriched KEGG pathway maps using Pathview package (v1.16.5) (http://bioconductor.org/packages/release/bioc/html/pathview.html) (64).
Reverse Transcription and Real-Time Quantitative PCR (RT-qPCR)
Reverse transcription was performed at 37°C for 2 h in a total volume of 10 μl reaction solution containing 3 μg RNA, 0.5 μg oligo (dT), 1X MMLV RT buffer, 0.5 mM each dNTP, 0.1 mM dithiothreitol, and 100 U M-MLV reverse transcriptase (Invitrogen, Carlsbad, CA). To validate the RNA-seq data, we determined the expression levels of some selected genes by RT-qPCR in both PG and PV follicles, including fshr, cyp19a1a, inha, inhbaa, notch3, amh, gadd45ga, lpl, march4, ube2ql1, zp2.1, and zp2.5. The expression levels were normalized to that of the housekeeping gene ef1a. The standard for each gene was prepared by PCR amplification of cDNA fragments with specific primers (Supplemental Table 1). After purification from the gel, the amplicons were quantitated by electrophoresis along with O'GeneRuler DNA Ladder Mix (Thermo Scientific, Waltham, MA), and the copy number of DNA molecule was calculated. The real-time qPCR assay was performed on the CFX96 Real-time PCR Detection System (Bio-Rad, Hercules, CA) and repeated twice.
Data Analysis
The mRNA level of each target gene was normalized to the internal control ef1a. All values were expressed as the mean ± SEM (n = 3), and the data were analyzed by t-test using Prism 6 on Macintosh OS X (GraphPad Software, San Diego, CA).
Results
Differential Gene Expression Profiles in PG and PV Follicles
The PG and PV follicles were isolated as described. The PG follicles were small and transparent with diameters less than 100 μm, whereas the PV follicles were much larger and less transparent. The PV follicles collected for analysis ranged from 180 to 250 μm in diameter. In histological sections, the PV follicles were bigger with numerous cortical alveoli in the oocytes and a continual layer of follicle cells surrounding the oocyte (Figure 1). Sequence analysis of mRNAs from PG and PV follicles showed that on average 84% reads were aligned to the zebrafish reference genome. The DEGs between PG and PV follicles are shown in four groups in the Volcano plot based on log2(fold change) and q-value (Benjamini-Hochberg adjusted p-value) [-log10(q-value)]: (1) significant for both fold change and q-value, (2) significant for fold change but nonsignificant for q-value, (3) nonsignificant for fold change but significant for q-value, and (4) nonsignificant for both fold change and q-value (Figure 2A). A large number of genes showed significant differential expression between the two stages (Figure 2B). To confirm the RNA-seq result, we evaluated 6 up-regulated genes (cyp19a1a, fshr, inha, inhbaa, notch3, and amh) and 6 down-regulated genes (gadd45ga, lpl, march4, ube2ql1, zp2.1, and zp2.5) (marked in Figure 2A) using RT-qPCR. The result showed that all these genes exhibited significant difference in their expression between PG and PV stages, fully agreeing with the RNA-seq result (Figure 3).
Figure 1. Morphology of zebrafish PG and PV follicles. (A) Ovarian fragments with PG and PV follicles. (B) Histological morphology of PG and PV follicles. (C) Isolated PG follicles (< 100 μm in diameter). (D) Isolated PV follicles (180 to 250 μm in diameter).
Figure 2. Volcano plot and heatmap for DEGs between PG and PV follicles. (A) Volcano plot for DGEs in PG/PV transition, showing distribution of significance [-log10(q-value)] vs. fold change [log2(fold change)] for all genes. The genes selected for qPCR validation are highlighted with black circles and gene symbols marked, including 6 down-regulated and 6 up-regulated genes. Genes in green are significant DEGs [q-value < = 0.05 and log2(fold change) > = 1]. (B) Heatmap of significant DGEs between PG and PV stage. Regularized log transformed (rlog) count matrix was generated using DeSeq2 package. Significant DEGs were used to plot the heatmap of rlog counts.
Figure 3. Validation of RNA-seq data by real-time qPCR on 12 significant DEGs highlighted in the Volcano plot, including 6 up-regulated and 6 down-regulated DEGs. (A) Up-regulated DGEs (PG vs. PV): cyp19a1a, amh, fshr, notch3, inhbaa, and inha. (B) Down-regulated DGEs (PG vs. PV): gadd45ga, lpl, march4, ube2ql1, zp2.1, and zp2.5. All the selected DGEs showed the same expression pattern in both RNA-seq and qPCR analyses. **P < 0.01; ***P < 0.001.
To classify the genes expressed in PG and PV follicles, including total expressed genes and DEGs (both up- and down-regulated genes), we generated the Venn diagram by an online software (http://bioinformatics.psb.ugent.be/webtools/Venn/). In total, 20,173 genes were detected with 17,074 expressed in both PG and PV follicles, and 1,313 and 1,786 detected in PG and PV only respectively. Among the genes expressed in both stages, 1,917 and 835 genes were up- and down-regulated in PV follicles respectively with statistical significance in both fold change and q-value. For 1,786 genes specifically expressed at PV stage, 110 showed significance in in both fold change and q-value whereas 24 PG-specific down-regulated genes showed significance (Figure 4A; for details see Supplemental Table 2). All these DEGs with significant change were subject to further GO analysis.
Figure 4. Venn diagram and scatter plots of enriched GO terms for up-regulated and down-regulated DEGs in PG/PV transition. (A) Venn diagram showing distribution of significant DEGs in different classes, including groups for expressed genes in PG and PV follicles, and up- and down-regulated gene groups. Among the genes expressed in both PG and PV follicles, 17,074 genes showed no significant difference in expression, and 1,917 and 835 were significantly up- and down-regulated, respectively, from PG to PV stage. A total of 1,786 genes were detected in PV follicles only, and 110 were significantly up-regulated in PG/PV transition. On the other hand, a total of 1,313 genes were detected in PG follicles only, and 24 were significantly down-regulated genes in PG/PV transition. (B) Scatter plot for enriched GO terms for up-regulated DEGs in PG/PV transition. (C) Scatter plot for enriched GO terms for down-regulated DEGs in PG/PV transition. Rich factor is the ratio of number of DEGs for particular GO term over the number of genes in background for that GO term. The size and color of the circles in scatter plots are scaled with respect to the number of DEGs and p-value respectively.
GO Analysis for DEGs in PG and PV Follicles
After analysis of gene expression profiles in PG and PV follicles, the DEGs with statistical significance were subject to further GO enrichment analysis, which covered all 2,027 (1917 + 110) up-regulated and 859 (835 + 24) down-regulated genes mentioned above. Among the three GO domains, we focused our analysis on the domain of Biological Process (BP).
The analysis revealed significant enrichment of a large number of GO terms during PG-PV transition. As expected, many GO terms associated with fundamental biological processes were enriched for up-regulated genes during follicle activation, including transcription (e.g., DNA-templated transcription), translation (e.g., amino acid transmembrane transport), metabolism (e.g., ATP synthesis coupled electron transport), and key cellular activities (cell growth, differentiation, migration, adhesion, and death). In addition, a variety of key signal transduction processes were enhanced, as indicated by the enrichment of up-regulated genes in GO terms associated with steroid hormone mediated signaling pathway, small GTPase mediated signal transduction, receptor tyrosine kinase signaling pathway, receptor serine/threonine kinase signaling pathway, TGF-β receptor signaling pathway, BMP signaling pathway, response to cAMP, response to cytokine, smoothened signaling pathway, ephrine receptor signaling pathway, and cellular response to hormone stimulus (Figure 4B). The up-regulated genes of some of these GO terms and their expression levels are illustrated in Figure 5. Interestingly, a series of GO terms associated with angiogenesis and vascularization were also up-regulated during follicle activation, including angiogenesis, sprouting angiogenesis, lymphangiogenesis, lymph vessel development, vasculogenesis, and blood vessel development (Figures 4B, 6). While a large number of GO terms were up-regulated during PG-PV transition, a number of GO terms were enriched for down-regulated genes, and they are mostly associated with chromatin structure (e.g., nucleosome assembly and chromatin silencing), egg-sperm interaction (e.g., egg coat formation, acrosome reaction, and sperm binding), RNA silencing (e.g., piRNA metabolic process and gene silencing by RNA), meiosis (e.g., meiotic cell cycle and meiotic DNA double-strand break formation), and proteolysis (Figure 4C). The representative GO terms enriched for down-regulated genes and the expression levels of these genes are plotted in Figure 7.
Figure 5. Graphic illustration of expression levels of various DEGs up-regulated in different biological processes related to cell signaling in PV follicles. The GO terms chosen for plotting were those with P < 0.1. *P < 0.05; **P < 0.01; ***P < 0.001.
Figure 6. Graphic illustration of expression levels of various DEGs up-regulated in different biological processes related to angiogenesis and blood vessel formation in PV follicles. The GO terms chosen for plotting were those with P < 0.1. *P < 0.05; **P < 0.01; ***P < 0.001.
Figure 7. Graphic illustration of expression levels of various DEGs down-regulated in different biological processes related to egg coat formation and egg-sperm interaction, piRNA metabolism and RNA silencing, and proteolysis in PV follicles. The GO terms chosen for plotting were those with P < 0.1. *P < 0.05; **P < 0.01; ***P < 0.001.
Pathway Analysis for DEGs in PG and PV Follicles
In addition to enrichment analysis for GO terms, we also performed KEGG pathway analysis to identify pathways potentially involved in follicle activation. Nineteen pathways were found to be significantly enriched for genes up-regulated during PG-PV transition (p < 0.05). Interestingly, most of these pathways are involved in signal transduction, including phosphatidylinositol signaling pathway, VEGF signaling pathway, hedgehog signaling pathway, FoxO signaling pathway, MAPK signaling pathway, TGF-β signaling pathway, Jak-STAT signaling pathway, ErbB signaling pathway, ECM-receptor interaction and cytokine-cytokine receptor interaction (Figure 8A). Several metabolic pathways were also enriched for up-regulated DEGs, including glycerophospholipid metabolism and glycosaminoglycan biosynthesis. Interestingly, a large number of up-regulated genes belonged to the pathway of endocytosis (Figure 8A). Compared with up-regulated pathways, fewer pathways were found to be down-regulated including tight junction and tryptophan metabolism (Figure 8B). The enriched KEGG pathways for both up- and down-regulated genes are shown in the Supplemental Figures 1, 2, respectively.
Figure 8. Scatter plots of enriched KEGG pathways. (A) Enriched KEGG pathways for up-regulated DEGs in PG/PV transition. (B) Enriched KEGG pathways for down-regulated DEGs in PG/PV transition. Rich factor is the ratio of number of DEGs for particular KEGG pathway over the number of genes in background for that KEGG pathway. The size and color of the circles in the plots are scaled with respect to the number of DEGs and p-value respectively.
Expression of Oocyte- and Follicle Cell-Specific Genes During Follicle Activation
Ovarian follicles consist of oocytes and somatic follicle cells. According to our published and unpublished data, the oocyte-specific genes often showed little changes in expression during follicle activation, whereas those expressed in the follicle cells are often up- or down-regulated. To provide further evidence for this observation, we examined the expression changes of genes specifically expressed in the oocytes and follicle cells. Since this issue has not been systematically studied in the zebrafish, the information about the localization of gene expression in the two compartments was based on previous studies in the zebrafish and a dataset available in the sheep (65), which was the first comprehensive RNA-seq study in mammals on transcriptomes of oocyte and granulosa cells in early follicles.
In agreement with our observations, most of genes expressed in the oocytes showed stable expression from PG to PV stage, including some well-known oocyte factors such as kita (Kit receptor A) and ybx1 (Y box binding protein 1). Some oocyte-specific genes showed decreased expression, including rxfp2b (Insl3 receptor), spo11 (initiator of meiotic double stranded breaks), and members of zona pellucida (ZP) family proteins, whereas some were increased such as abcb4 (a member of the ATP-binding cassette transporter family) and slc7a3a/b (cationic amino acid transporter 3). In contrast, a large number of genes expressed in the follicle cells displayed significant changes in their expression during the PG-PV transition and most were up-regulated. Many of the top up-regulated genes are related to steroidogenesis, synthesis of estrogens and their actions, and signaling of various endocrine and paracrine factors, including nr5a1a/b (steroidogenic factor 1), cyp19a1a (ovarian aromatase) and its transcription stimulator foxl2 (forkhead box L2), esr1 and esr2b (estrogen receptors), fshr (FSH receptor), notch3 (neurogenic locus notch homolog protein 3), inha and inhbaa (inhibin and activin subunits), acvrl1 (a type I receptor of TGF-β superfamily), wnt4a (a ligand of frizzled family receptors), egfra and erbb2 (receptors of EGF family), bmpr2a/b (type II receptors of BMPs), and kitlga (Kit ligand A) (Figure 9).
Figure 9. Heatmaps of gene expression change during PG/PV transition with respect to cell type specificity in the ovarian follicle, viz. oocytes and follicle cells. The genes in red indicates those that increased expression from PG to PV whereas those in blue decreased expression in the process.
Discussion
The follicle development in the zebrafish consists of two major phases, i.e., PG and SG. The SG phase is further divided into previtellogenic (PV), early vitellogenic (EV), mid-vitellogenic (MV), late vitellogenic (LV), and full-grown (FG) stages according to follicle size and morphological characteristics (3, 66, 67). Although zebrafish is now being commonly used as a model for studying fish reproduction and its folliculogenesis in the ovary is well characterized at morphological level, the molecular mechanisms underlying follicle activation or PG-PV transition remain largely unknown. Using NGS-based RNA-seq, the present study demonstrated, for the first time, a comprehensive gene expression profile during ovarian follicle activation in the zebrafish. Our data showed that although most of the genes detected were common between the two follicle stages, many genes showed stage specificity in expression. Our analysis revealed a large number of candidate genes that were potentially involved in the activation of ovarian follicles or PG-PV transition, which also marks puberty onset when it occurs for the first time in zebrafish life cycle.
The quality of RNA-seq analysis was validated by real-time qPCR on twelve selected genes with six up- and six down-regulated from PG to PV. These genes were chosen based on either their positions in the Volcano plot or previous studies. All twelve genes showed consistent results in RNA-seq and qPCR analyses. The up-regulation of cyp19a1a, fshr, amh, inha and inhbaa is consistent with previous studies (11, 23, 36, 68) and the increased expression of notch3 agrees with our unpublished data (Lau and Ge, unpublished). A recent RNA-seq study in the coho salmon also demonstrated increased expression of cyp19a1, fshr, amh and inha in 11-KT-induced PG-PV transition (19). Interestingly, the expression of zona pellucida (Zp/zp) genes (zp2 and zp3) showed significant reduction during PG-PV transition. This seems contradictory to the role of Zp proteins in follicle development because the chorion, which consists of Zp proteins, is one of the hallmarks for follicle activation or PG-PV transition (23). This pattern of zp expression, however, agrees well with previous studies in both zebrafish and other species. The transcripts of zp2 and zp3 were abundantly expressed in the PG oocytes, and their levels decreased significantly while entering the PV stage. This was also confirmed by in situ hybridization, transgenesis and promoter analysis (69–71). Furthermore, their expression was independent of estrogen signaling (70). This has been confirmed in other fish species such as coho salmon (57).
In addition to the aforementioned genes, further analysis of the RNA-seq data identified a large number of biological processes and pathways that are potentially involved in controlling follicle activation, including those concerning transcription, translation, metabolism, and various cellular activities such as proliferation, differentiation, migration, and death. These are expected because follicle activation involves significant structural remodeling. Herein we will focus our discussion on the following aspects: signal transduction, angiogenesis, endocytosis, glycosaminoglycan biosynthesis, RNA silencing, and proteolysis.
Signal Transduction
Follicle activation is a critical step of early folliculogenesis, which involves dramatic changes in both structure and function (1, 72). A large body of evidence has accumulated in recent years in both mammals and fish that follicle activation or recruitment involves extensive regulation by endocrine and paracrine factors (20, 73, 74), which activate a variety of signaling pathways (75). As master hormones that control gonadal development and function, the pituitary gonadotropins (FSH and LH) signal mostly through cAMP/PKA pathway (76, 77). Loss of FSH receptor (fshr) in the zebrafish led to complete arrest of folliculogenesis at PG stage (4), suggesting an important role for gonadotropin signaling at follicle activation. Our demonstration of increased fshr expression (67) and enrichment of biological process response to cAMP during PG-PV transition provides supportive evidence for such regulation.
Together with gonadotropins, steroid hormones are believed to be essential for normal folliculogenesis, in particular the estrogen signaling pathway (78, 79). The involvement of gonadal steroids in follicle activation is supported in the present study by enriched GO terms of steroid hormone mediated signaling pathway and cellular response to estrogen stimulus. Among the genes significantly up-regulated were nuclear estrogen receptors (esr1 and esr2b), androgen receptor (ar) and progestin receptor (pgr). The increased expression of nuclear estrogen receptors agrees well with our recent study in the zebrafish (80) and the reports in the coho salmon that E2 had potent stimulatory effect on the growth of PV follicles and accumulation of cortical alveoli (18, 19). In mammals, androgen signaling has been reported to be essential for normal folliculogenesis (81), and the disruption of androgen signaling resulted in defective folliculogenesis (82). In fish, lines of evidence from different species also point to the importance of androgens in early follicle development, especially growth of PG follicles and PG-PV transition (14, 16–19). Using RNA-seq approach, a recent study in the coho salmon demonstrated significant changes in ovarian transcriptome after exposure to 11-KT. A variety of genes were up-regulated by 11-KT, including genes associated with steroidogenesis and steroid action (19), similar to our data in this study. In Japanese huchen and common carp, both progestin (DHP, 17α, 20β-dihydroxy-4-pregnen-3-one) and E2 have been shown to regulate early oogenesis (83). A recent study in the zebrafish showed that the disruption of pgr led to female infertility due to failed ovulation (84). It would be interesting to examine if the loss of pgr has any impact on early follicle development. Taken together, steroid-mediated signaling pathway might be one of the pivotal mechanisms that regulate PG-PV transition in the zebrafish ovary.
In addition to gonadotropins and steroids, a large number of local ovarian paracrine factors especially peptide growth factors have been reported to play roles in controlling mammalian folliculogenesis especially during early stages (76, 85), including families of IGF (86), EGF (87), TGF-β (88), activin (89), BMP including BMP15 and GDF9 (90), Kit (91), and Notch (92). These growth factors and their signaling components also exist in the zebrafish ovary (22, 24–27, 30, 34, 93, 94). The importance of these ovarian paracrine factors in controlling zebrafish follicle activation is strongly supported by our RNA-seq data. A series of GO terms concerning signal transduction pathways, mostly those in growth factor signaling, were enriched in up-regulated genes, including transmembrane receptor protein tyrosine kinase signaling pathway, TGF-β receptor signaling pathway, transmembrane receptor protein serine/threonine kinase signaling pathway, BMP signaling pathway, and Notch signaling pathway. Consistent with GO enrichment analysis, the KEGG pathway analysis also suggested seminal importance for various cellular signaling pathways in follicle activation. Most pathways with significant enrichment for up-regulated genes turned out to be signal transduction pathways, including phosphatidylinositol, VEGF, TGF-β, hedgehog, FoxO, MAPK, Jak-STAT, and ErbB signaling pathways. Most of these pathways have been implicated in ovarian follicle activation or recruitment in mammals (95–97). Among these pathways, the phosphatidylinositol or PI3K/PTEN/Akt signaling system is well known to be particularly important for follicle recruitment or activation in mammals (98). The defects in this signaling pathway are closely associated with female infertility (99, 100). In summary, our transcriptomics data provide strong support to traditional views that ovarian follicle activation or recruitment involves pituitary gonadotropins, steroid hormones and local growth factors.
Angiogenesis
In addition to signal transduction pathways as described above, one of the most interesting discoveries of the present study was that a large number of biological processes related to angiogenesis or blood vessel formation were enriched in up-regulated genes during PG-PV transition, including angiogenesis, sprouting angiogenesis, vasculogenesis, blood vessel development, blood vessel morphogenesis, lymph vessel formation, response to hypoxia, and regulation of VEGF receptor signaling pathway. The VEGF signaling pathway was also significantly enriched in KEGG pathway analysis as described above. This result suggests that follicle activation or PG-PV transition must involve active angiogenesis or formation of blood vessels, which is likely to ensure sufficient supply of nutrients, hormones and oxygen to support fast growing follicles.
In mammals, the ovary is one of the few organs where physiological angiogenesis occurs in adult, and ovarian angiogenesis has been well documented to be essential for normal follicle growth and maturation (101, 102). A recent study in infantile rats demonstrated a significant increase in VEGF production during preantral-antral transition, further supporting a role for angiogenesis during follicle activation (103). In fish, the oocyte growth involves incorporation of hepatic vitellogenin from the maternal blood stream, making angiogenesis and sufficient blood supply a critical condition for normal folliculogenesis. So far there have been very few studies on ovarian angiogenesis in fish models. In tilapia, both VEGF and its receptor VEGFR2 (Flk-1) were shown to be expressed in the granulosa cells with increasing levels during follicle growth, whereas angiopoietin-2 (Ang-2) and its receptor Tie-2 were expressed in the theca and granulosa cells respectively with expression mostly observed in early SG follicles (104). This agrees well with our transcriptome data in the present study.
Endocytosis
Different from follicle growth in mammals, which mostly involves proliferation of somatic granulosa and theca cells surrounding the oocyte, the growth of follicles in fish and other oviparous species, especially in the SG stage, is largely due to phenomenal growth of the oocytes because of accumulation of hepatic yolk protein vitellogenin (105). The uptake of vitellogenin from the blood stream by the growing oocytes is mediated by receptor-dependent endocytosis (106, 107). In agreement with this, our pathway analysis revealed endocytosis pathway as one of the major pathways that were significantly enriched with a large number of up-regulated genes. This strongly suggests that although the exogenous vitellogenesis does not occur at PV stage, the increased expression of endocytosis-related genes prepares the follicles for vitellogenin uptake in the subsequent phase of oocyte growth.
Glycosaminoglycan Biosynthesis
The glycosaminoglycans (GAGs) are the most abundant heteropolysaccharides in the body and they are also important components of proteoglycans. The GAG biosynthesis includes two pathways for heparan sulfate/heparin and chondroitin sulfate/dermatan sulfate synthesis respectively. Interestingly, the genes involved in these pathways were significantly increased in PV follicles. It was recently reported that ovarian GAGs promoted gonadotropin-induced follicle development in the mouse, and these GAGs might act, at least partially, by stimulating angiogenesis (108), which agrees with the increased angiogenic activities discussed above. Earlier studies demonstrated that the synthesis of proteoglycans with core proteins and GAGs changed significantly during ovarian folliculogenesis regulated by FSH (109). These studies and our data all suggest important roles for GAGs during folliculogenesis, especially follicle activation.
Although the exact roles of GAGs in PG-PV transition remain unclear, we speculate that they may be related to the biogenesis of the cortical alveoli (CA) or cortical vesicles/granules in the oocytes, which is a major hallmark of PV stage as compared to PG (66, 110). The cortical alveoli play an important role in cortical reaction during fertilization (111), and their contents contain rich glycosylated materials (112) including glycoproteins such as polysialoglycoprotein (113, 114) and GAG-like substances (115).
RNA Silencing
Our GO and KEGG analyses showed that most enriched biological processes and pathways contained up-regulated genes during PG-PV transition, suggesting increased biological activities in the process. However, some GO terms were significantly enriched for down-regulated genes, including piRNA metabolic process and gene silencing by RNA. piRNAs (Piwi-interacting RNAs) are a special group of small RNAs (24–31 nt long) that are associated with the Piwi proteins, a subgroup of the Argonaute family. These small RNAs and Piwi are specifically expressed in the germ cells to target transposons and other genes (116). Evidence has accumulated that piRNAs play important roles in germ cell development, embryogenesis and genome stability (117) although the importance has not been supported by Piwi knockout study in the mouse (118).
In addition to piRNA, microRNA-mediated RNA silencing has been reported to play important roles in ovarian folliculogenesis (119). In a recent transcriptome study, we characterized the profiles of microRNAs during the PG-PV transition in the zebrafish ovary. We identified a 13-microRNA expression signature that showed significant changes during the transition (120). Although the role of non-coding small RNAs in zebrafish folliculogenesis remains largely unknown, the enrichment of both piRNA metabolic process and gene silencing by RNA in down-regulated genes during the PG-PV transition strongly suggests an important role for RNA interference or silencing in follicle activation. This would be an interesting direction for research in the future.
Proteolysis
Although it was not as significant as other GO terms, the proteolysis process was enriched by a large number of down-regulated genes during PG-PV transition. This is interesting and has promoted us to hypothesize that an active protein turnover may occur in PG and/or early PV follicles; however, the activity of protein degradation decreases when the follicles approach the end of PV stage (our sampling stage) prior to entering vitellogenic growth. Our unpublished proteomics analysis also showed a dramatic change in protein profile during PG-PV transition in contrast to the relatively consistent and stable profiles afterwards in vitellogenic growth (Lau and Ge, unpublished data).
Oocytes vs. Follicle Cells
In addition to biological process and pathway enrichment analyses during PG-PV transition, we also looked at the expression changes of genes that are expressed in the two different compartments of ovarian follicles, i.e., oocytes and follicle cells. Since the information regarding expression localization within the follicles is limited in the zebrafish, our analysis was mostly based on previously published data in the zebrafish and a dataset available in the sheep with the assumption that most genes also have similar spatial distribution in the zebrafish follicle (65). Interestingly, as we observed in our previous studies, most oocyte-expressed genes were stable in their expression during follicle activation except a few; in contrast, a large number of genes in the follicle cells showed dramatic changes in expression, mostly up-regulated, during PG-PV transition. This suggests that the regulation of oocyte-expressed genes may likely occur at post-transcriptional level whereas those in the follicle cells are mostly controlled at the transcription level. The observation that many up-regulated genes that showed the most significant changes in the follicle cells were related to signaling of endocrine and paracrine factors strongly suggests that the somatic follicle cells are the site where various signaling molecules act and interact to control folliculogenesis.
In summary, the transcriptome analysis in the present study for differentially expressed genes during follicle activation provides critical insights into the candidate molecules involved in early folliculogenesis, and in particular, the regulatory networks that play roles in triggering and driving follicle activation and puberty onset. As expected, a large number of signal transduction pathways were activated during PG-PV transition or follicle activation, which may be responsible for signaling of endocrine hormones and paracrine factors such as gonadotropins, steroids, and growth factors. In addition, a variety of pathways concerning angiogenesis or blood vessel formation were also activated during follicle activation, indicating increased demand for supplying nutrients, hormones, and oxygen to the growing follicles. Other pathways responsible for specific events or processes such as biogenesis of cortical alveoli and vitellogenin uptake were also significantly enriched. The follicle activation or PG-PV transition may also involve other regulatory mechanisms such as RNA silencing and proteolysis. With the powerful genome editing technologies (TALEN and CRISPR/Cas9) available now in the zebrafish model, the functions of novel genes identified can be further validated by targeted loss-of-function approach.
Ethics Statement
This study was carried out in accordance with the recommendations by the Legislative Council of Macao Special Administrative Region under Article 71 (1) of the Basic Law. The protocol was approved by the Research Ethics Committee of the University of Macau.
Author Contributions
BZ designed and executed all the experiments, writing the manuscript. LP: data analysis. YC: data analysis. WG: experimental design, writing the manuscript.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
This study was supported by grants from the University of Macau (MYRG2014-00062-FHS, MYRG2015-00227-FHS, MYRG2016-00072-FHS, and CPG2014-00014-FHS) and The Macau Fund for Development of Science and Technology (FDCT114/2013/A3 and FDCT/089/2014/A2) to WG. We thank the Genomics and Bioinformatics Core facility of the Faculty of Health Sciences for experimental and technical support on sequencing and data analysis. All the bioinformatics analyses were performed on High Performance Computing (HPC) facility supported by Information and Communication Technology Office (ICTO) of the University of Macau. We thank Jacky Chan and William Pang for their technical support on HPC.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2018.00593/full#supplementary-material
References
1. Mcgee EA, Hsueh AJ. Initial and cyclic recruitment of ovarian follicles. Endocr Rev. (2000) 21:200–14. doi: 10.1210/edrv.21.2.0394
2. Kezele PR, Ague JM, Nilsson E, Skinner MK. Alterations in the ovarian transcriptome during primordial follicle assembly and development. Biol Reprod. (2005) 72:241–55. doi: 10.1095/biolreprod.104.032060
3. Ge W. Paracrine control of fish ovarian follicle development and function. In: Garcia-Ayala A, Penalver JM, Chaves-Pozo E, editors. Recent Advances in Fish Reproduction Biology (Research Signpost) (2010) p. 141–73.
4. Zhang Z, Lau SW, Zhang L, Ge W. Disruption of zebrafish follicle-stimulating hormone receptor (fshr) but not luteinizing hormone receptor (lhcgr) gene by TALEN leads to failed follicle activation in females followed by sexual reversal to males. Endocrinology (2015) 156:3747–62. doi: 10.1210/en.2015-1039
5. Taranger GL, Carrillo M, Schulz RW, Fontaine P, Zanuy S, Felip A, et al. Control of puberty in farmed fish. Gen Comp Endocrinol. (2010) 165:483–515. doi: 10.1016/j.ygcen.2009.05.004
6. Chen W, Ge W. Ontogenic expression profiles of gonadotropins (fshb and lhb) and growth hormone (gh) during sexual differentiation and puberty onset in female zebrafish. Biol Reprod. (2012) 86:73. doi: 10.1095/biolreprod.111.094730
7. Chen W, Ge W. Gonad differentiation and puberty onset in the zebrafish: evidence for the dependence of puberty onset on body growth but not age in females. Mol Reprod Dev. (2013) 80:384–92. doi: 10.1002/mrd.22172
8. Khoo KH. The histochemistry and endocrine control of vitellogenesis in goldfish ovaries. Can J Zool. (1979) 57:617–26. doi: 10.1139/z79-073
9. Swanson P. Salmon gonadotropins: reconciling old and new ideas. In: Scott AP, Sumpter JP, Kime DE, Rolfe MS, editors. Proceedings of the Fourth International Symposium on the Reproductive Physiology of Fish (Sheffield: FishSymp91) (1991), p. 2–7.
10. Campbell B, Dickey J, Beckman B, Young G, Pierce A, Fukada H, et al. Previtellogenic oocyte growth in salmon: relationships among body growth, plasma insulin-like growth factor-1, estradiol-17beta, follicle-stimulating hormone and expression of ovarian genes for insulin-like growth factors, steroidogenic-acute regulatory protein and receptors for gonadotropins, growth hormone, and somatolactin. Biol Reprod. (2006) 75:34–44. doi: 10.1095/biolreprod.105.049494
11. Kwok HF, So WK, Wang Y, Ge W. Zebrafish gonadotropins and their receptors: I. Cloning and characterization of zebrafish follicle-stimulating hormone and luteinizing hormone receptors - evidence for their distinct functions in follicle development. Biol Reprod. (2005) 72:1370–81. doi: 10.1095/biolreprod.104.038190
12. So WK, Kwok HF, Ge W. Zebrafish gonadotropins and their receptors: II. Cloning and characterization of zebrafish follicle-stimulating hormone and luteinizing hormone subunits - their spatial-temporal expression patterns and receptor specificity. Biol Reprod. (2005) 72:1382–96. doi: 10.1095/biolreprod.104.038216
13. Zhang Z, Zhu B, Ge W. Genetic analysis of zebrafish gonadotropin (FSH and LH) functions by TALEN-mediated gene disruption. Mol Endocrinol. (2015) 29:76–98. doi: 10.1210/me.2014-1256
14. Lokman PM, George KA, Divers SL, Algie M, Young G. 11-Ketotestosterone and IGF-I increase the size of previtellogenic oocytes from shortfinned eel, Anguilla australis, in vitro. Reproduction (2007) 133:955–67. doi: 10.1530/REP-06-0229
15. Rohr DH, Lokman PM, Davie PS, Young G. 11-Ketotestosterone induces silvering-related changes in immature female short-finned eels, Anguilla australis. Comp Biochem Physiol A Mol Integr Physiol. (2001) 130:701–14. doi: 10.1016/S1095-6433(01)00402-0
16. Kortner TM, Rocha E, Arukwe A. Androgenic modulation of early growth of Atlantic cod (Gadus morhua L.) previtellogenic oocytes and zona radiata-related genes. J Toxicol Environ Health A. (2009). 72:184–95. doi: 10.1080/15287390802539020
17. Kortner TM, Rocha E, Arukwe A. Previtellogenic oocyte growth and transcriptional changes of steroidogenic enzyme genes in immature female Atlantic cod (Gadus morhua L.) after exposure to the androgens 11-ketotestosterone and testosterone. Comp Biochem Physiol A Mol Integr Physiol. (2009) 152:304–13. doi: 10.1016/j.cbpa.2008.11.001
18. Forsgren KL, Young G. Stage-specific effects of androgens and estradiol-17β on the development of late primary and early secondary ovarian follicles of coho salmon (Oncorhynchus kisutch) in vitro. Biol Reprod. (2012) 87:64. doi: 10.1095/biolreprod.111.098772
19. Monson C, Forsgren K, Goetz G, Harding L, Swanson P, Young G. A teleost androgen promotes development of primary ovarian follicles in coho salmon and rapidly alters the ovarian transcriptome. Biol Reprod. (2017) 97:731–45. doi: 10.1093/biolre/iox124
20. Ge W. Intrafollicular paracrine communication in the zebrafish ovary: the state of the art of an emerging model for the study of vertebrate folliculogenesis. Mol Cell Endocrinol. (2005) 237:1–10. doi: 10.1016/j.mce.2005.03.012
21. Spiliotis BE. Growth hormone insufficiency and its impact on ovarian function. Ann NY Acad Sci. (2003) 997:77–84. doi: 10.1196/annals.1290.009
22. Wu T, Patel H, Mukai S, Melino C, Garg R, Ni X, et al. Activin, inhibin, and follistatin in zebrafish ovary: expression and role in oocyte maturation. Biol Reprod. (2000) 62:1585–92. doi: 10.1095/biolreprod62.6.1585
23. Wang Y, Ge W. Developmental profiles of activin βA, βB, and follistatin expression in the zebrafish ovary: evidence for their differential roles during sexual maturation and ovulatory cycle. Biol Reprod. (2004) 71:2056–64. doi: 10.1095/biolreprod.104.032649
24. Kohli G, Hu S, Clelland E, Di Muccio T, Rothenstein J, Peng C. Cloning of transforming growth factor-β1 (TGF-β1) and its type II receptor from zebrafish ovary and role of TGF-β1 in oocyte maturation. Endocrinology (2003) 144:1931–41. doi: 10.1210/en.2002-0126
25. Clelland E, Kohli G, Campbell RK, Sharma S, Shimasaki S, Peng C. Bone morphogenetic protein-15 in the zebrafish ovary: complementary deoxyribonucleic acid cloning, genomic organization, tissue distribution, and role in oocyte maturation. Endocrinology (2006) 147:201–9. doi: 10.1210/en.2005-1017
26. Liu L, Ge W. Growth differentiation factor 9 and its spatiotemporal expression and regulation in the zebrafish ovary. Biol Reprod. (2007) 76:294–302. doi: 10.1095/biolreprod.106.054668
27. Li CW, Ge W. Spatiotemporal expression of bone morphogenetic protein family ligands and receptors in the zebrafish ovary: a potential paracrine-signaling mechanism for oocyte-follicle cell communication. Biol Reprod. (2011) 85:977–86. doi: 10.1095/biolreprod.111.092239
28. Dranow DB, Hu K, Bird AM, Lawry ST, Adams MT, Sanchez A, et al. Bmp15 Is an oocyte-produced signal required for maintenance of the adult female sexual phenotype in zebrafish. PLoS Genet. (2016) 12:e1006323. doi: 10.1371/journal.pgen.1006323
29. Wang Y, Ge W. Cloning of epidermal growth factor (EGF) and EGF receptor from the zebrafish ovary: evidence for EGF as a potential paracrine factor from the oocyte to regulate activin/follistatin system in the follicle cells. Biol Reprod. (2004) 71:749–60. doi: 10.1095/biolreprod.104.028399
30. Tse AC, Ge W. Spatial localization of EGF family ligands and receptors in the zebrafish ovarian follicle and their expression profiles during folliculogenesis. Gen Comp Endocrinol. (2010) 167:397–407. doi: 10.1016/j.ygcen.2009.09.012
31. Chung CK, Ge W. Epidermal growth factor differentially regulates activin subunits in the zebrafish ovarian follicle cells via diverse signaling pathways. Mol Cell Endocrinol. (2012) 361:133–42. doi: 10.1016/j.mce.2012.03.022
32. Li J, Liu Z, Wang D, Cheng CH. Insulin-like growth factor 3 is involved in oocyte maturation in zebrafish. Biol Reprod. (2011) 84:476–86. doi: 10.1095/biolreprod.110.086363
33. Li J, Chu L, Sun X, Liu Y, Cheng CH. IGFs mediate the action of LH on oocyte maturation in zebrafish. Mol Endocrinol. (2015) 29:373–83. doi: 10.1210/me.2014-1218
34. Zhou R, Yu SM, Ge W. Expression and functional characterization of intrafollicular GH-IGF system in the zebrafish ovary. Gen Comp Endocrinol. (2016) 232:32–42. doi: 10.1016/j.ygcen.2015.11.018
35. Yao K, Ge W. Kit system in the zebrafish ovary: evidence for functional divergence of two isoforms of Kit (Kita and Kitb) and Kit ligand (Kitlga and Kitlgb) during folliculogenesis. Biol Reprod. (2010) 82:1216–26. doi: 10.1095/biolreprod.109.082644
36. Rodriguez-Mari A, Yan YL, Bremiller RA, Wilson C, Canestro C, Postlethwait JH. Characterization and expression pattern of zebrafish anti-Mullerian hormone (amh) relative to sox9a, sox9b, and cyp19a1a, during gonad development. Gene Expr Patterns (2005) 5:655–67. doi: 10.1016/j.modgep.2005.02.008
37. Qin M, Zhang Z, Song W, Wong QW, Chen W, Shirgaonkar N, et al. Roles of Figla/figla in juvenile ovary development and follicle formation during zebrafish gonadogenesis. Endocrinology (2018). doi: 10.1210/en.2018-00648. [Epub ahead of print].
38. Lau ES, Zhang Z, Qin M, Ge W. Knockout of zebrafish ovarian aromatase gene (cyp19a1a) by TALEN and CRISPR/Cas9 leads to all-male offspring due to failed ovarian differentiation. Sci Rep. (2016) 6:37357. doi: 10.1038/srep37357
39. Yin Y, Tang H, Liu Y, Chen Y, Li G, Liu X, et al. Targeted disruption of aromatase reveals dual functions of cyp19a1a during sex differentiation in zebrafish. Endocrinology (2017) 158:3030–41. doi: 10.1210/en.2016-1865
40. Nakamoto M, Shibata Y, Ohno K, Usami T, Kamei Y, Taniguchi Y, et al. Ovarian aromatase loss-of-function mutant medaka undergo ovary degeneration and partial female-to-male sex reversal after puberty. Mol Cell Endocrinol. (2018) 460:104–22. doi: 10.1016/j.mce.2017.07.013
41. Li MH, Yang HH, Li MR, Sun YL, Jiang XL, Xie QP, et al. Antagonistic roles of Dmrt1 and Foxl2 in sex differentiation via estrogen production in tilapia as demonstrated by TALENs. Endocrinology (2013) 154:4814–25. doi: 10.1210/en.2013-1451
42. Yang YJ, Wang Y, Li Z, Zhou L, Gui JF. Sequential, divergent, and cooperative requirements of Foxl2a and Foxl2b in ovary development and maintenance of zebrafish. Genetics (2017) 205:1551–72. doi: 10.1534/genetics.116.199133
43. Xie QP, He X, Sui YN, Chen LL, Sun LN, Wang DS. Haploinsufficiency of SF-1 causes female to male sex reversal in Nile tilapia, Oreochromis niloticus. Endocrinology (2016) 157:2500–14. doi: 10.1210/en.2015-2049
44. Nishimura T, Sato T, Yamamoto Y, Watakabe I, Ohkawa Y, Suyama M, et al. foxl3 is a germ cell-intrinsic factor involved in sperm-egg fate decision in medaka. Science (2015) 349:328–31. doi: 10.1126/science.aaa2657
45. Chu L, Li J, Liu Y, Cheng CH. Gonadotropin signaling in zebrafish ovary and testis development: insights from gene knockout study. Mol Endocrinol. (2015) 29:1743–58. doi: 10.1210/me.2015-1126
46. Kocabas AM, Crosby J, Ross PJ, Otu HH, Beyhan Z, Can H, et al. The transcriptome of human oocytes. Proc Natl Acad Sci USA. (2006) 103:14027–32. doi: 10.1073/pnas.0603227103
47. Koks S, Velthut A, Sarapik A, Altmae S, Reinmaa E, Schalkwyk LC, et al. The differential transcriptome and ontology profiles of floating and Cumulus granulosa cells in stimulated human antral follicles. Mol Hum Reprod. (2010) 16:229–40. doi: 10.1093/molehr/gap103
48. Labrecque R, Fournier E, Sirard MA. Transcriptome analysis of bovine oocytes from distinct follicle sizes: Insights from correlation network analysis. Mol Reprod Dev. (2016) 83:558–69. doi: 10.1002/mrd.22651
49. Bonnet A, Servin B, Mulsant P, Mandon-Pepin B. Spatio-temporal gene expression profiling during in vivo early ovarian folliculogenesis: integrated transcriptomic study and molecular signature of early follicular growth. PLoS ONE (2015) 10:e0141482. doi: 10.1371/journal.pone.0141482
50. Magalhaes-Padilha DM, Geisler-Lee J, Wischral A, Gastal MO, Fonseca GR, Eloy YR, et al. Gene expression during early folliculogenesis in goats using microarray analysis. Biol Reprod. (2013) 89:19. doi: 10.1095/biolreprod.112.106096
51. Tosser-Klopp G, Le Cao KA, Bonnet A, Gobert N, Hatey F, Robert-Granie C, et al. A pilot study on transcriptome data analysis of folliculogenesis in pigs. Animal (2009) 3:393–401. doi: 10.1017/S1751731108003479
52. Xu X, Zhao X, Lu L, Duan X, Qin H, Du X, et al. Transcriptomic analysis of different stages of pigeon ovaries by RNA-sequencing. Mol Reprod Dev. (2016) 83:640–8. doi: 10.1002/mrd.22670
53. Tingaud-Sequeira A, Chauvigne F, Lozano J, Agulleiro MJ, Asensio E, Cerda J. New insights into molecular pathways associated with flatfish ovarian development and atresia revealed by transcriptional analysis. BMC Genomics (2009) 10:434. doi: 10.1186/1471-2164-10-434
54. Von Schalburg KR, Rise ML, Brown GD, Davidson WS, Koop BF. A comprehensive survey of the genes involved in maturation and development of the rainbow trout ovary. Biol Reprod. (2005) 72:687–99. doi: 10.1095/biolreprod.104.034967
55. Jia Y, Nan P, Zhang W, Wang F, Zhang R, Liang T, et al. Transcriptome analysis of three critical periods of ovarian development in Yellow River carp (Cyprinus carpio). Theriogenology (2018) 105:15–26. doi: 10.1016/j.theriogenology.2017.08.027
56. Knoll-Gellida A, Andre M, Gattegno T, Forgue J, Admon A, Babin PJ. Molecular phenotype of zebrafish ovarian follicle by serial analysis of gene expression and proteomic profiling, and comparison with the transcriptomes of other animals. BMC Genomics (2006) 7:46. doi: 10.1186/1471-2164-7-46
57. Luckenbach JA, Iliev DB, Goetz FW, Swanson P. Identification of differentially expressed ovarian genes during primary and early secondary oocyte growth in coho salmon, Oncorhynchus kisutch. Reprod Biol Endocrinol. (2008) 6:2. doi: 10.1186/1477-7827-6-2
58. Wang Y, Ge W. Gonadotropin regulation of follistatin expression in the cultured ovarian follicle cells of zebrafish, Danio rerio. Gen Comp Endocrinol. (2003) 134:308–15. doi: 10.1016/S0016-6480(03)00275-2
59. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. (2013) 14:R36. doi: 10.1186/gb-2013-14-4-r36
60. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122
61. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. (2016) 11:1650–67. doi: 10.1038/nprot.2016.095
62. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
63. Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, et al. The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. (2007) 8:R183. doi: 10.1186/gb-2007-8-9-r183
64. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics (2013) 29:1830–1. doi: 10.1093/bioinformatics/btt285
65. Bonnet A, Cabau C, Bouchez O, Sarry J, Marsaud N, Foissac S, et al. An overview of gene expression dynamics during early ovarian folliculogenesis: specificity of follicular compartments and bi-directional dialog. BMC Genomics (2013) 14:904. doi: 10.1186/1471-2164-14-904
66. Selman K, Wallace RA, Sarka A, Qi XP. Stages of oocyte development in the zebrafish, Brachydanio rerio. J Morphol. (1993) 218:203–24. doi: 10.1002/jmor.1052180209
67. Zhou R, Tsang AH, Lau SW, Ge W. Pituitary adenylate cyclase-activating polypeptide (PACAP) and its receptors in the zebrafish ovary: evidence for potentially dual roles of PACAP in controlling final oocyte maturation. Biol Reprod. (2011) 85:615–25. doi: 10.1095/biolreprod.111.091884
68. Poon SK, So WK, Yu X, Liu L, Ge W. Characterization of inhibin alpha subunit (inha) in the zebrafish: evidence for a potential feedback loop between the pituitary and ovary. Reproduction (2009) 138:709–19. doi: 10.1530/REP-09-0198
69. Onichtchouk D, Aduroja K, Belting HG, Gnugge L, Driever W. Transgene driving GFP expression from the promoter of the zona pellucida gene zpc is expressed in oocytes and provides an early marker for gonad differentiation in zebrafish. Dev Dyn. (2003) 228:393–404. doi: 10.1002/dvdy.10392
70. Liu X, Wang H, Gong Z. Tandem-repeated zebrafish zp3 genes possess oocyte-specific promoters and are insensitive to estrogen induction. Biol Reprod. (2006) 74:1016–25. doi: 10.1095/biolreprod.105.049403
71. Mold DE, Dinitz AE, Sambandan DR. Regulation of zebrafish zona pellucida gene activity in developing cocytes. Biol Reprod. (2009) 81:101–10. doi: 10.1095/biolreprod.108.071720
72. Mcgee EA, Raj RS. Regulators of ovarian preantral follicle development. Semin Reprod Med. (2015) 33:179–84. doi: 10.1055/s-0035-1552584
73. Hsueh AJ, Kawamura K, Cheng Y, Fauser BC. Intraovarian control of early folliculogenesis. Endocr Rev. (2015) 36:1–24. doi: 10.1210/er.2014-1020
74. Atwood CS, Vadakkadath Meethal S. The spatiotemporal hormonal orchestration of human folliculogenesis, early embryogenesis and blastocyst implantation. Mol Cell Endocrinol. (2016) 430:33–48. doi: 10.1016/j.mce.2016.03.039
76. Erickson GF, Shimasaki S. The physiology of folliculogenesis: the role of novel growth factors. Fertil Steril. (2001) 76:943–9. doi: 10.1016/S0015-0282(01)02859-X
77. Ge W. Gonadotropins and their paracrine signaling network in the zebrafish ovary. Fish Physiol Biochem. (2005) 31:209–14. doi: 10.1007/s10695-006-0026-1
78. Drummond AE. The role of steroids in follicular growth. Reprod Biol Endocrinol. (2006) 4:16. doi: 10.1186/1477-7827-4-16
79. Prossnitz ER, Arterburn JB, Smith HO, Oprea TI, Sklar LA, Hathaway HJ. Estrogen signaling through the transmembrane G protein-coupled receptor GPR30. Annu Rev Physiol. (2008) 70:165–90. doi: 10.1146/annurev.physiol.70.113006.100518
80. Liu KC, Lau SW, Ge W. Spatiotemporal expression analysis of nuclear estrogen receptors in the zebrafish ovary and their regulation in vitro by endocrine hormones and paracrine factors. Gen Comp Endocrinol. (2017) 246:218–25. doi: 10.1016/j.ygcen.2016.12.011
81. Kimura S, Matsumoto T, Matsuyama R, Shiina H, Sato T, Takeyama K, et al. Androgen receptor function in folliculogenesis and its clinical implication in premature ovarian failure. Trends Endocrinol Metab. (2007) 18:183–9. doi: 10.1016/j.tem.2007.04.002
82. Hu YC, Wang PH, Yeh SY, Wang RS, Xie C, Xu QQ, et al. Subfertility and defective folliculogenesis in female mice lacking androgen receptor. Proc Natl Acad Sci USA. (2004) 101:11209–14. doi: 10.1073/pnas.0404372101
83. Miura C, Higashino T, Miura T. A progestin and an estrogen regulate early stages of oogenesis in fish. Biol Reprod. (2007) 77:822–8. doi: 10.1095/biolreprod.107.061408
84. Zhu Y, Liu D, Shaner ZC, Chen S, Hong W, Stellwag EJ. Nuclear progestin receptor (Pgr) knockouts in zebrafish demonstrate role for pgr in ovulation but not in rapid non-genomic steroid mediated meiosis resumption. Front Endocrinol. (2015) 6:37. doi: 10.3389/fendo.2015.00037
85. Monget P, Monniaux D. Growth factors and the control of folliculogenesis. J Reprod Fertil. (1995) 49:321–33.
86. Monget P, Bondy C. Importance of the IGF system in early folliculogenesis. Mol Cell Endocrinol. (2000) 163:89–93. doi: 10.1016/S0303-7207(99)00244-0
87. Park JY, Su YQ, Ariga M, Law E, Jin SL, Conti M. EGF-like growth factors as mediators of LH action in the ovulatory follicle. Science (2004) 303:682–4. doi: 10.1126/science.1092463
88. Drummond AE. TGFβ signalling in the development of ovarian function. Cell Tissue Res. (2005) 322:107–15. doi: 10.1007/s00441-005-1153-1
89. Findlay JK, Drummond AE, Dyson M, Baillie AJ, Robertson DM, Ethier JF. Production and actions of inhibin and activin during folliculogenesis in the rat. Mol Cell Endocrinol. (2001) 180:139–44. doi: 10.1016/S0303-7207(01)00521-4
90. Paradis F, Novak S, Murdoch GK, Dyck MK, Dixon WT, Foxcroft GR. Temporal regulation of BMP2, BMP6, BMP15, GDF9, BMPR1A, BMPR1B, BMPR2 and TGFBR1 mRNA expression in the oocyte, granulosa and theca cells of developing preovulatory follicles in the pig. Reproduction (2009) 138:115–29. doi: 10.1530/REP-08-0538
91. Hutt KJ, Mclaughlin EA, Holland MK. Kit ligand and c-Kit have diverse roles during mammalian oogenesis and folliculogenesis. Mol Hum Reprod. (2006) 12:61–9. doi: 10.1093/molehr/gal010
92. Vanorny DA, Mayo KE. The role of Notch signaling in the mammalian ovary. Reproduction (2017) 153:R187–204. doi: 10.1530/REP-16-0689
93. Yao K, Ge W. Spatial distribution and receptor specificity of zebrafish Kit system - evidence for a Kit-mediated bi-directional communication system in the preovulatory ovarian follicle. PLoS ONE (2013) 8:e56192. doi: 10.1371/journal.pone.0056192
94. Chen W, Liu L, Ge W. Expression analysis of growth differentiation factor 9 (Gdf9/gdf9), anti-mullerian hormone (Amh/amh) and aromatase (Cyp19a1a/cyp19a1a) during gonadal differentiation of the zebrafish, Danio rerio. Biol Reprod. (2017) 96:401–13. doi: 10.1095/biolreprod.116.144964
95. Knight PG, Glister C. TGF-β superfamily members and ovarian follicle development. Reproduction (2006) 132:191–206. doi: 10.1530/rep.1.01074
96. Du XY, Huang J, Xu LQ, Tang DF, Wu L, Zhang LX, et al. The proto-oncogene c-src is involved in primordial follicle activation through the PI3K, PKC and MAPK signaling pathways. Reprod Biol Endocrinol. (2012) 10:58. doi: 10.1186/1477-7827-10-58
97. Zhao Y, Zhang Y, Li J, Zheng N, Xu X, Yang J, et al. MAPK3/1 participates in the activation of primordial follicles through mTORC1-KITL signaling. J Cell Physiol. (2018) 233:226–37. doi: 10.1002/jcp.25868
98. Reddy P, Liu L, Adhikari D, Jagarlamudi K, Rajareddy S, Shen Y, et al. Oocyte-specific deletion of Pten causes premature activation of the primordial follicle pool. Science (2008) 319:611–3. doi: 10.1126/science.1152257
99. Hu MW, Wang ZB, Schatten H, Sun QY. New understandings on folliculogenesis/oogenesis regulation in mouse as revealed by conditional knockout. J Genet Genomics (2012) 39:61–8. doi: 10.1016/j.jgg.2012.01.004
100. Makker A, Goel MM, Mahdi AA. PI3K/PTEN/Akt and TSC/mTOR signaling pathways, ovarian dysfunction, and infertility: an update. J Mol Endocrinol. (2014) 53:R103–18. doi: 10.1530/JME-14-0220
102. Augustin HG. Vascular morphogenesis in the ovary. Baillieres Best Pract Res Clin Obstet Gynaecol. (2000) 14:867–82. doi: 10.1053/beog.2000.0132
103. Torres-Ortiz MC, Gutierrez-Ospina G, Gomez-Chavarin M, Murcia C, Alonso-Morales RA, Perera-Marin G. The presence of VEGF and Notch2 during preantral-antral follicular transition in infantile rats: anatomical evidence and its implications. Gen Comp Endocrinol. (2017) 249:82–92. doi: 10.1016/j.ygcen.2017.05.006
104. Santos FC, Silva JF, Boeloni JN, Teixeira E, Turra EM, Serakides R, et al. Morphological and immunohistochemical characterization of angiogenic and apoptotic factors and the expression of thyroid receptors in the ovary of tilapia Oreochromis niloticus in captivity. Pesq Vet Bras. (2015) 35:371–6. doi: 10.1590/S0100-736X2015000400010
105. Hiramatsu N, Todo T, Sullivan CV, Schilling J, Reading BJ, Matsubara T, et al. Ovarian yolk formation in fishes: molecular mechanisms underlying formation of lipid droplets and vitellogenin-derived yolk proteins. Gen Comp Endocrinol. (2015) 221:9–15. doi: 10.1016/j.ygcen.2015.01.025
106. Patino R, Sullivan CV. Ovarian follicle growth, maturation, and ovulation in teleost fish. Fish Physiol Biochem. (2002) 26:57–70. doi: 10.1023/A:1023311613987
107. Mizuta H, Luo W, Ito Y, Mushirobira Y, Todo T, Hara A, et al. Ovarian expression and localization of a vitellogenin receptor with eight ligand binding repeats in the cutthroat trout (Oncorhynchus clarki). Comp Biochem Physiol B Biochem Mol Biol. (2013) 166:81–90. doi: 10.1016/j.cbpb.2013.07.005
108. Sato E. Intraovarian control of selective follicular growth and induction of oocyte maturation in mammals. Proc Jpn Acad Ser B Phys Biol Sci. (2015) 91:76–91. doi: 10.2183/pjab.91.76
109. Salustri A, Camaioni A, Di Giacomo M, Fulop C, Hascall VC. Hyaluronan and proteoglycans in ovarian follicles. Hum Reprod Update (1999) 5:293–301.
110. Wallace RA, Selman K. Cellular and dynamic aspects of oocyte growth in teleosts. Amer Zool. (1981) 21:325–43.
111. Sun QY. Cellular and molecular mechanisms leading to cortical reaction and polyspermy block in mammalian eggs. Microsc Res Tech. (2003) 61:342–8. doi: 10.1002/jemt.10347
112. Hoodbhoy T, Talbot P. Mammalian cortical granules: contents, fate, and function. Mol Reprod Dev. (1994) 39:439–48.
113. Kitajima K, Inoue Y, Inoue S. Polysialoglycoproteins of Salmonidae fish eggs - complete structure of 200-kDa polysialoglycoprotein from the unfertilized eggs of rainbow trout (Salmo gairdneri). J Biol Chem. (1986) 261:5262–9.
114. Selman K, Wallace RA, Barr V. Oogenesis in Fundulus heteroclitus .V. The relationship of yolk vesicles and cortical alveoli. J Exp Zool. (1988) 246:42–56.
115. Wessel GM, Brooks JM, Green E, Haley S, Voronina E, Wong J, et al. The biology of cortical granules. Int Rev Cytol. (2001) 209:117–206. doi: 10.1016/S0074-7696(01)09012-X
116. Iwasaki YW, Siomi MC, Siomi H. PIWI-interacting RNA: its biogenesis and functions. Annu Rev Biochem. (2015) 84:405–33. doi: 10.1146/annurev-biochem-060614-034258
117. Watanabe T, Lin H. Posttranscriptional regulation of gene expression by Piwi proteins and piRNAs. Mol Cell (2014) 56:18–27. doi: 10.1016/j.molcel.2014.09.012
118. Rosenkranz D, Han CT, Roovers EF, Zischler H, Ketting RF. Piwi proteins and piRNAs in mammalian oocytes and early embryos: from sample to sequence. Genom Data (2015) 5:309–13. doi: 10.1016/j.gdata.2015.06.026
119. Fitzgerald JB, George J, Christenson LK. Non-coding RNA in ovarian development and disease. Adv Exp Med Biol. (2016) 886:79–93. doi: 10.1007/978-94-017-7417-8_5
Keywords: transcriptomics, RNA-seq, follicles, folliculogenesis, puberty, zebrafish
Citation: Zhu B, Pardeshi L, Chen Y and Ge W (2018) Transcriptomic Analysis for Differentially Expressed Genes in Ovarian Follicle Activation in the Zebrafish. Front. Endocrinol. 9:593. doi: 10.3389/fendo.2018.00593
Received: 26 June 2018; Accepted: 18 September 2018;
Published: 11 October 2018.
Edited by:
Yong Zhu, East Carolina University, United StatesReviewed by:
Chun Peng, York University, CanadaGraham Young, University of Washington, United States
Balasubramanian Senthilkumaran, University of Hyderabad, India
Copyright © 2018 Zhu, Pardeshi, Chen and Ge. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence:Wei Ge, weige@umac.mo; gezebrafish@gmail.com