Androgens stimulate primary ovarian development in Vertebrate. Japanese eels underwent operation to sample the pre- and post-treated ovarian tissues from the same individual. Ovarian phenotypic or genotypic data were mined in a pair. A correlation between the initial ovarian status (determined by kernel density estimation (KDE), presented as a probability density of oocyte size) and the consequence of androgen (17MT) treatment (change in ovary) has been showed. The initial ovarian status appeared to be important to influence ovarian androgenic sensitivity. The initial ovary was important to the outcomes of androgen treatments, and ePAV (expression presence-absence variation) is existing in Japanese eel by analyze DEGs; core, unique, or accessory genes were identified, the sensitivities of initial ovaries were correlated with their gene expression profiles. We speculated the importance of genetic differential expression on the variations of phenotypes by 17MT, and transcriptomic approach seems to allow extracting multiple layers of genomic data.
- Anguilla japonica
- kernel density estimation
- primary ovary
There are different responsive abilities of different eels on piscine pituitary extracts (PPE) in artificial maturation manipulation. This phenomenon has been commented: “the maturity status of parental females at the beginning of artificial maturation trials is also an important factor affecting the success of maturation experiments”  as well as suggested that silver eels with more developed oocytes (mean of oocyte diameter is 0.262 mm) required a shorter time to produce mature eggs than did cultivated eels (mean of that is 0.177 mm) , the maturity of morphological changes (silvering, as a pubertal metamorphosis in eels) might have a close correlation with ovarian development (oocyte growth). The physiological status of the eels just before PPE injection has been postulated to be a reason for this difference . In addition, in most circumstances during artificial eel maturation, unsynchronized development stages of ovary are a problem, androgens have been demonstrated to synchronize oocyte development in the eel ( references therein). While, in practice, the stimulatory effect of androgens (or PPE) on eel ovarian development was often not uniform or not conserved among different individuals, and the stimulatory effects might depend on the initial status of ovarian development [3, 4]. Artificial maturation in eels appears to be paralogue to ‘personalized medicine’ in human diseases. Identifying a set of gene biomarkers to predict or estimate the outcomes of eel artificial maturation is interesting and is demanded.
2. From phenotypes
Indeed, the developmental stages of oocytes in Japanese eels were classified by the mean of oocyte diameter as follows: oil droplet stage for follicles <200 μm; early vitellogenic stage for follicles 200–400 μm; mid-vitellogenic stage for follicles 400–600 μm; late vitellogenic stage for follicles 600–800 μm; and migratory nucleus stage for follicles 800–1000 μm [4, 6].
Androgens stimulate the development of the preantral follicle in mammals . In fish, androgens participate in the promotion of primary and/or early secondary ovarian follicle growth (reviewed by ), the primary and/or early secondary ovarian follicle stages in the fish is analogous to the preantral follicle in mammals. The stimulatory effect of androgens has been reported in coho salmon (
Indeed, artificial maturation is often unsuccessful when the eel broodstocks come from aquafarms or have a larger body size . Whether, the hypotheses that the outcomes of artificial eel maturation are based on “the maturity status of parental females at the beginning of artificial maturation trials”  and “the physiological status of the eels just before hormone”  are accepted. Therefore, identifying sex and evaluating gonadal maturity are necessary.
A method to track ovarian development before- and after-androgen (17α-methyltestosterone, 17MT) treatment in the same eel has been developed . It must be emphasized that the experimental data came from the same individual; this traceability was the advantage, but the bias caused by a single sample could not be avoided. The pretreated ovary was used as a morphological or molecular baseline to compare the posttreatment ovary. Since the experimental eels and PPE were both from wild populations and subject to important individual variation, both of them are not consisted. This means the results based on two important variations if the different eels as well as PPE were employed simultaneously. To increase sample size can reduce the variation (to approach normal distribution), it can be speculated that the normal distribution is composed of many eels, the number of eels will be numerous, when one eel was selected, since time cannot be reversal, the normal distribution will collapse into one sample, then the destinations or outcomes of the selected one might be an another distribution (not always a normal distribution), but, apparently, there should be a closer correlation between the inputs and the outcomes. The data followed this protocol can be accumulated then be clustered, different patterns or clusters related with certain biological outcomes may be distinguished.
Eels were surgically operated to identify sex and to sample ovarian tissues, the operated eels recovered and the sutures were shed after 4–5 weeks [13, 5]. The manipulation did not suspend ovarian development under the subsequent artificial maturation regimen.
Primary ovarian sensitivity to androgens (17MT) is well documented; this characteristic was used to study the correlation between basal ovarian status and the outcome in Japanese eel. The lard grease- and Vaseline-based sustained-release mixtures were used
In the pretreatment eels, there was no significant (p > 0.05) difference among them on mean of oocyte diameter, the developmental stage of our experimental eels was the oil droplet stage (primary oocyte). After 4 weeks of treatment (posttreatment), The mean diameter of oocyte significantly (p < 0.05) increased in the same individual, and the mean of fold change was x1.25. Based on KDE and using 0.1 mm as the arbitrary boundary, the percentage of the oocyte diameter over 0.1 mm increased from 44–70% .
Three eels, with different fold change on oocyte diameter x1.30, x2.08, and x2.35 tagged with no 2003, no 1440, and no 2613, respectively, were chosen to present the different results of the17MT treatment, although the growth of all oocytes was significantly (p < 0.05) stimulated by LVMT, as revealed by pre- vs. posttreatment comparisons (Figure 1). Moreover, analyzed by KDE, and 0.1 mm was used as an arbitrary boundary to evaluate ovarian status, the percentage of the oocytes over 0.1 mm increased from 22–60% in eel no 2003, from 42–90% in eel no 1440, and from 40–89% in eel no 2613. It seemed that oocyte growth (development) was more synchronized in eel no 2613 than in eel no 1440 because the bandwidth in eel no 2613 (bandwidth = 0.009722) was smaller than that in eel no 1440 (bandwidth = 0.02096), which indicates a more dispersed trend in ovarian development in eel no 1440; these data corresponded to the standard deviations (SDs) of 0.0574 and 0.0352 for eel no 1440 and no 2613, respectively (Figure 2) .
KDE may be a way to mine the “maturity status” of the ovary. From plotted figure, if we chose 0.1 mm as the boundary, we can speculate that if the main peak was located on or right-hand to the boundary, better oocyte development results could be detected under the subsequent treatment condition with 17MT.
More than two peaks can identified in certain plotted figures, a shoulder instead of a smooth bell curve was also observed, hidden peak was likely present . Our observed data implied a wave-like on the development of Japanese eel oocytes ; coincidently, the Japanese eel might spawn more than once during a spawning season . Moreover, there is also an argument that if the Japanese eel is a not semelparous (monocyclic) fish .
The pretreated ovary was used as a morphological or molecular baseline to compare the posttreatment ovary in the same individual. The comparison was based on the pattern from a probability density function (KDE); the patterns represented ovarian status. Our results showed that a different ovarian status seemed to lead to a different artificial maturation results .
3. To genotypes
In Japanese eel, the correlations between the initial ovarian status, represented by the probability distributions of oocyte diameter, and the effects of 17MT treatment is shown . We continued to investigate the relationships between the initial ovarian gene expression profiles and the stimulatory effects of androgens on oocyte development in eels. As the activation and development of ovaries are sophisticated (reviewed by [17, 18]), myriad factors and ample genes interact with each other to control folliculogenesis . Following this concept, we studied this phenomenon based on the ovarian expression profiles of the initial ovaries using transcriptomic tools, high-throughput methods for the large-scale analysis of ovarian gene functions can be used to approach this question.
Gene set analysis (GSA) has the advantage of incorporating existing biological knowledge into the expression analysis (reviewed by , references therein). Gene Ontology (GO) terms can be used to define gene sets, and GO analysis is a standard way to define the relationships among expressed genes, thus enabling identification through the use of GSA. Both GO and COG (Cluster of Orthologous Genes) provide specific information about genes or gene products.
In teleosts, there is a lack of transcriptome analysis on early folliculogenesis (reviewed by ). In coho salmon (
The number of total trinity transcripts was 541,020, and 452,778 transcripts were derived from raw assembled transcripts but clustered with 95% similarity. The number of total trinity ‘genes’ was 336,973, while a total of 25,340 unigenes were distributed into the 3 categories with 37.56% in biological processes, 49.15% in molecular functions, and 13.28% in cellular components. Meanwhile, the top 10 genes in biological processes, cellular component, and molecular function were indicated. According to the annotation of COG (clusters of orthologous groups of proteins), these genes were classified into 25 different functional classes. The cluster for general function prediction (19,089, 24.77%) represented the largest group, followed by cell cycle control, cell division, chromosome partitioning (11,670, 15.14%), replication, recombination and repair (6,704; 8.70%), and unknown functions (5,286; 6.86%) .
The post-treatment ovarian tissue was compared with the initial ovarian tissue in the same eel. The results showed that certain pathways were significantly (
REViGO (http://revigo.irb.hr) was used to summarize lists of GO terms by finding a representative subset of the terms. The analysis was based on the post-treatment ovarian tissue compared with the initial ovary (e.g., 2003MT/2003). The total GO enrichment score of eels 2003, 2613, and 1440 was 184, 253, and 230, respectively. These results corresponded with the phenotypes. The outcomes from the GO enrichment analysis and phenotypes were mutually supported and appeared to have a positive correlation. The ovarian gene expression profiles from the same eel before-and after-treatment are investigated. The profiles of the heatmap or cluster indicated that eel 2613 and eel 1440 were under the same clade (cluster), and eel 2003 appeared to be isolated from both. And, the gene expression profiles were more discrete in eel 2613 between the initial ovary and post-treated ovary compared with those of the two other eels (Figure 3). Plotting the TMM expression of all samples and calculating the Pearson correlation coefficient, the coefficient was 0.34, 0.17, and 0.15 for eel 2003, 1440, and 2613, respectively, between the initial ovary and post-treated one in the same individual .
The Principal Component Analysis (PCA), the distance between any given two points represents their divergences in the chart, results are also shown: the distance between eels 1440MT and 2613MT was short on PC1 (X-axis), while eel 2003MT was isolated on the other side; eel 26130MT was closer to 2003MT while 1440MT was located from the other two on PC2 (Y-axis). The order of the intra-group distance was: 2613 > 1440 > 2003; the major shift was along PC1, and there was slight shift along PC2 except for eel 2613. Moreover, the results from the heatmap also indicated that eels 2613MT and 1440MT under the same clade (cluster), but eel 2003MT was isolated from both. The results are supported mutually (Figure 4) .
We preferred to mine the differences in gene expression profiles, thoroughly as these differences might provide the bases to set up biomarkers. To maximize the differences significantly, the cutoff condition was set with the TPM of the gene as zero in the initial ovarian tissues. We attempted to filter the differences between the sensitive group (2613MT and 1440MT) and insensitive one (2003MT), the TPM was set to be equal or greater than 0.01 (TPM ≥ 0.01) to amplify the differences; there were 467 genes exclusively in 2613MT and 518 genes exclusively in 1440MT. As there were many isoforms, the number of unique genes was 266 and 303 for 2613MT and 1440MT, respectively. These unique genes are involved in various biological pathways. There were many GO enrichment biological pathways in addition to housekeeping metabolic pathways, the neurogenesis pathways (axon guidance, eye morphogenesis, and axon development in 2613MT; sensory organ development, regulation of synapse organization, and regulation of nervous system development in 1440MT) appeared to overlap in these two eels. These data implied plausible roles of neurogenesis in the primary development of the ovary .
The cutoff levels were TPM ≥ 0.01 and TPM = 0 for the sensitive group and insensitive group, respectively, tried to filter data based on DEGs to maximize the differences. Our assumption was that the significant (≥ 0.01
Interestingly, several pathways implicated in neurogenesis or neuro-activity have been mined. We changed the cutoff-levels; curiously, the pathways implicating eyes or sensory organ development were found. The retina is the most accessible part of the vertebrate central nervous system, and the neural retina constitutes an excellent system with which to analyze key aspects of neurogenesis . Moreover, it is well recognized that the silver eel (the prepubertal stage) has larger eyes, and a similar situation might occur in the developing gonads. Our cutoff-based results led us to postulate that neurogenesis and/or neuro-activity is involved in primary ovarian development in Japanese eel. This hypothesis might pertain to the fact that in eel, even entering the silver stage, the reproductive endocrine system is blocked, but a slight development of ovary is observed.
Few reports have discussed the involvement of neurogenesis in primary ovarian development in fish. Meanwhile, in mammals, the development of ovarian nerves precedes the onset of folliculogenesis , and the neural activity might be an important factor in the regulation of follicular development before the ovary acquires responsiveness to gonadotropins , the developing ovary is regulated by direct sympathetic inputs which function in addition but complementarily to gonadotropin . On the other hand, the density of innervation may contribute to the selection of follicles for further development (reviewed by ). Evidence supports complementary control performed by sympathetic nerves innervating the ovary [29, 30]. Moreover, our non-
Indeed, our non-
Biomarkers are biological characteristics that predict treatment responses. The sensitivity of oocytes to 17MT is dependent on the ovarian status (not the mean size of the oocyte) . From our results, it seems not ideal to select certain genes based on DEGs to serve as biomarkers, instead of independent genes with significant DEGs, a set of genes in certain pathways to represent the ovarian status seems chosen. The initial ovarian status is fundamental to influencing how genotypes are translated into phenotypes by treatments. It is already useful to inspect the size of sampled oocytes, while our report indicated that, the gene expression pattern can ameliorate and improve the method based on phenotype for the eel artificial maturation. Our results also implied that the neurogenesis or neural activity in the initial ovary might be an important factor influencing the outcomes of hormonal treatments. Meanwhile, more experiments and data are required to support and verify this hypothesis.
4. From ePAV to phenotypes
The term “pan-genome” is the entirety of a species gene repertoire, it is defined as the complete gene set across strains rather than in a single individual. The concept of pan-genome is built upon the observation that genes often display presence-absence variations (PAV) . Pan-genome comprises of the “core” genes share by all strains within a given species, “dispensable” or “accessory” genes are only present in some strains, “unique” genes are strain-specific . PAV is an important source of genetic divergence and diversity  and has been suggested to contribute to phenotypic variation of agronomic importance in various crops . With the completion of reference genome sequencing, our knowledge of genomic variation increased, it is apparent that a single reference sequence is insufficient to represent the extent of genomic variation, resulting in the introduction and growth of the pan-genome concept . In nowadays, PAV has been commonly observed across the Tree of Life .
The pan-transcriptome can be defined by recalling the concept of the pan-genome. It reflects the set of all the RNA present in a specie or in a single organism, transcripts present in every individual of a taxonomic group are called core genes, while transcripts absent in at least a single individual are called accessory (dispensable) genes . Indeed, a larger part of the genes are expression presence-absence variation (ePAV) or are even totally absent (genomic PAV; gPAV) . PAV at the genomic level would be reflected in the transcriptome ePAV. The ePAV not only reflect genomic structural variation, but also the variations in genetic and epigenetic regulatory elements .
As a correlation between the initial ovarian status and the phenotypic outcomes of 17MT treatment have been demonstrated , and the genetic background of initial ovary seems important . We ameliorated our results that the ovarian transcriptome was re-mapped against the Japanese genomic scaffolds, a pair as the unit to re-analyze DEGs. Certain genes were absent in one pair but present in the others, the presence genes were also up- or down-regulated, the genes of various sets were analyzed for KEGG pathway enrichment. Together up, we speculated the existence of presence-absence variant (PAV) or expression presence-absence variation (ePAV) in Japanese eel.
The 22199 pan-ovarian transcripts were derived from ovarian tissues of 12 eels with different maturation stages. In a previous report, the number of total trinity “genes” was 336,973 , while there were 22,199 pan-ovarian transcripts after re-mapped against the Japanese genomic scaffolds. This significant discrepancy shows the importance of specific genomic scaffolds on transcriptomics. The number of functional genes was countered under the criterion that a gene was detected in at least one of six samples (three pairs), and the gene was eliminated when this gene was not detected in any of the six samples. The number of core genes was 15466; there were 980, 501, and 1097 accessory genes for the intersection between 1440s and 2003s, between 2003s and 2613 s, and between 2613 s and 1440s, respectively. The numbers of unique genes were 1335, 796, and 724 for the 1440s, 2003s, and 2613 s, respectively (Figure 5). we emphasize PAVs here. We cannot confirm that the PAVs were stage-specific or inducible and would be detectable (expressed) in the later stages, since the ovaries of this stage remained far from maturation.
Regarding the phenotype, there was no significant difference in mean ovarian diameter among the initial ovarian tissues (Figure 1) , but all genotype data from both PCA and heatmaps indicate a closer distance between 1440MT (the initial ovarian tissue) and 2613MT, while 2003MT was more isolated from others (Figures 3 and 4) , these observations imply that there were genetic variations. Accordingly, we speculate that the expression of the absent genes may be up-regulated (from zero to detectable) by the treatment or in the following stage; if so, these absent genes are not presence-absence variations. It has been demonstrated that, not all absent genes were developmental- or stage-specific, which indirectly supports the concept of ePAV in Japanese eel .
It is essential to characterize the ePAV genes and their possible functions , the core genes has been speculated to appear universally over-represented by house-keeping functions and essential metabolic processes for organisms, and accessory genes are often associated with communication, virulence, and defense responses (, references therein] or adaptive functions in plants (e.g., stress responses ). It has also been demonstrated that accessory genes exhibit higher rates of polymorphism than core genes [44, 45]. We studied possible functions of the unique genes of PAVs in each pair by KEGG pathway enrichment testing, since phenotypes are determined or controlled by constellations of pathways. The list of unique presence genes of each pair was submitted to Molas (http://molas.iis.sinica.edu.tw/jpeel/) to test for the KEGG pathway enrichment, and the pathways were marked when the enrichment p-value <0.05; these pathways were clustered into seven categories in the KEGG pathway database. Indeed, only approximately 36% of unique genes were hit in the KEGG pathways (the hit rates were 37%, 36%, and 36% for 1440s, 2003s and 2613 s, respectively). These results indicate that the functions of transcripts in eels remained largely unknown and non-annotated, while the results also indicated that 1440s was clearly more similar to 2003.
All signaling pathways belong to the category of “Environmental Information Processing” in KEGG, and the numbers of pathways of each pair in this category were 12, 7, and 5 for 1440s, 2003s and 2613, respectively. There were 3 pathways at the intersection of those 3 pairs: Rap1 signaling pathway, cGMP-PKG signaling pathway, and cAMP signaling pathway. This result implies the importance of these signaling pathways in basal ovarian development, and the importance of differential pathways in ovarian development was shown. Moreover, we believe that the differentially expressed genes are integrated into certain pathways to determine or control phenotypes. It has been demonstrated that the functions of genes related to neuronal activities or neurogenesis appears important in the outcomes of androgen treatments . This postulation is supported the existing of ePAV in Japanese eel.
Our results, from Japanese eel, implied that the initial ovary status might be an important factor influencing the outcomes of androgen treatment, and the genes related with neuronal activities or neurogenesis seemed to play an essential role in the positive effect; certain genes seemed absent in one individual but present in the others by speculating that there is a presence-absence variant (PAV) or expression presence-absence variation (ePAV) in Japanese eel, a primitive fish diverged from other bony fishes in the basal lineage of the Teleost. Most pathways involved by ePAV were belong the endocrine system and nervous system. These results signify the importance of genetic differential expression on the variations of phenotypes by androgen, and a transcriptomic approach appears to enable extracting multiple layers of genomic data. More elegant experimental designs fit the biological model, and more transcriptomic data are required to address this question in the future.
We would like to thank Mr. Yung-Lin Huang and Mr. Yi-Sen Lai, their encouragements were appreciated sincerely. Mr. Guan-Ru Chen (Freshwater Aquaculture Research Center, Fisheries Research Institute, C.O.A., Taiwan) is also appreciated for his help on the manipulation. This research was partly funded by the Ministry of Science and Technology (MOST107-2321-B002-057), Taiwan.
Conflict of interest
The authors declare that they have no conflict of interest.
This study was partly funded by the Ministry of Science and Technology (MOST107–2321-B002–057), Taiwan.
All software application employed are online open access.