Transcriptome Profiling of the Ovarian Cells at the Single-Cell Resolution in Adult Asian Seabass

Single-cell RNA sequencing (scRNA-seq) is widely adopted for identifying the signature molecular markers or regulators in cells, as this would benefit defining or isolating various types of cells. Likewise, the signature transcriptome profile analysis at the single cell level would well illustrate the key regulators or networks involved in gametogenesis and gonad development in animals; however, there is limited scRNA-seq analysis on gonadal cells in lower vertebrates, especially in the sexual reversal fish species. In this study, we analyzed the molecular signature of several distinct cell populations of Asian seabass adult ovaries through scRNA-seq. We identified five cell types and also successfully validated some specific genes of germ cells and granulosa cells. Likewise, we found some key pathways involved in ovarian development that may concert germline-somatic interactions. Moreover, we compared the transcriptomic profiles across fruit fly, mammals, and fish, and thus uncovered the conservation and divergence in molecular mechanisms that might drive ovarian development. Our results provide a basis for studying the crucial features of germ cells and somatic cells, which will benefit the understandings of the molecular mechanisms behind gametogenesis and gonad development in fish.


INTRODUCTION
Sexually reproducing organisms transmitted their genetic information from one generation to the next via gametogenesis, a physiological procedure initiated with the primordial germ-cell (PGC) specification and migration (Kunwar and Lehmann, 2003). Traditionally, light and electron microscopy have been used to discriminate PGCs from somatic cells during early embryonic stages and examine their morphological characteristics (Xu et al., 2010). In recent years, based on molecular and cellular techniques, several identified germ cell markers such as vasa, dazl, and ziwi are known to play essential roles in germ cell development across the animal kingdom (Kehkooi et al., 2009;Leu and Draper, 2010;Aduma et al., 2019). However, the identification of germ cells in many species was ambiguous due to the absence of germ cells or somatic marker genes, especially for female germ cells in ovary.
The ovary is a highly conserved transcriptionally tissue across animal taxa from invertebrates to vertebrates (Edwards, 1965;Zagalsky et al., 1967;Pedersen and Peters, 1968;Thomson et al., 2010). Female germ cell development is a complex multifactorial regulated process in that embryonic germ cells are driven to become oogonia stem cells through a series of statuses including cell proliferation, apoptosis, and differentiation, then subsequently generate ootid through oogenesis (Sánchez and Smitz, 2012). Moreover, the female germ cell fate mainly relies on the ovarian environments established by somatic cells rather than the sex chromosomes in germ cells (Evans et al., 1977). Likewise, oogenesis has been extensively documented via the anatomical and histological analysis (Grier et al., 2017;Teh et al., 2018;Pierson, 2019); however, the molecular basis of germ cells is still largely unclear in fish species. Therefore, it is necessary to precisely define the various cell types in order to well understand the mechanism underlying ovarian development.
Single cell RNA sequencing (scRNA-seq) approaches have recently been applied to study the embryo development and gametogenesis due to its high accuracy in detecting the variations across cell lineages and heterogeneity or predict developmental trajectories (Li L. et al., 2017;Wagner et al., 2018). In contrast to the traditional transcriptome analysis, transcriptome profiling at the single-cell level crosses the limitation of obscured molecular heterogeneity among and between varied gametogenic cells, which provides comprehensive understanding and significant new biological insights (Suzuki et al., 2019). Developmental trajectories during zebrafish embryogenesis via scRNA-seq analysis provides a new approach to reconstruct the specification trajectories of many developmental systems without the need for prior knowledge of gene expression profiles, fatemaps, or lineage trajectories (Farrell et al., 2018). In human adult ovaries, the maps of the molecular signature of growing and regressing follicular populations constructed using scRNA-seq highlight the significance of mapping whole adult organs at the single-cell level, and this is also paramount to understand the mechanisms of (in)fertility . scRNA-seq analysis in human oocytes and corresponding somatic follicles provided new insights into understandings on the crucial features of transcriptional machinery, transcription factor networks, and reciprocal interactions between oocytes and somatic follicular compartments (Zhang et al., 2018). In adult Drosophila testis, the transcriptome profiles of genetic novelties in various cell types across spermatogenesis was elucidated by scRNA-seq analysis, offering significant clues for understanding how the testis maintains its core reproductive function while being a hotbed of evolutionary innovation (Witt et al., 2019). The scRNA-seq analysis has been widely applied to studying the germ cells or gonadal development in organisms in recent years. However, up to now, most investigations on ovarian transcriptome focus mainly on limited species, including human (Liu et al., 2009;Zhang et al., 2018), monkey , mouse (Dean, 2002;Zhao et al., 2020), zebrafish , and fly (Slaidina et al., 2020), which revealed the remarkable interspecies variations and conservations during ovarian development across phyla.
Lates calcarifer, also named Asian sea bass or barramundi, is an important economic species owing to its delicious taste and high nutritional value and has been successfully cultivated in saltwater throughout the Asian-Pacific region for more than 20 years (Mohd-Yusof et al., 2010;Talpur, 2014). Asian sea bass is hermaphroditic with male-first-matured that matures into a functional male at 2-4 years and then reverses into a female at subsequent spawning season. However, little is known about its molecular and cellular mechanisms behind the sex reverse of male-to-female in animals, including fish. In this study, scRNA-seq analysis was performed and aimed to profile the gene expression of all cell types in Asian seabass ovary. Here, five distinct types of cells including one germ cell and four somatic cells were classified. Likewise, a series of new candidate marker genes were detected, and this would facilitate identifying and isolating the various types of ovarian cells in Asian sea bass in our future studies. Also, several potential biological processes and signaling pathways involved in germ cells development and differentiation were figured out. Additionally, a comparative analysis on gene expression profiles during folliculogenesis were conducted, and some conserved genes and species-specific genes associated with the development of female germ cells were identified between mammals and fish. Here, we identified both germ cells and ovarian somatic cells, and successfully demonstrated the first single-cell transcriptomic atlas of seabass adult ovary. Therefore, this study would pave the way for the further extensive investigations on gene expression profiles in fish ovarian cells, and shed new insights into understanding the conservation and divergence in the mechanisms behind the germ cells differentiation across phyla.

Single-Cell Isolation of Ovary
Adult Asian sea bass at age of 4 + years old were collected from a fish farm at Hainan Haikou. The ovarian cells of various types were isolated for single cell transcriptomics analysis as previously described (Xu et al., 2018;Hochane et al., 2019;Witt et al., 2019). Briefly, after being washed 2-3 times with 1 × phosphate buffer saline (PBS) and cut into pieces, the ovarian tissue was successively digested with 1 mg/ml Collagenase Type I (LS004196, Worthington) and 0.25% Trypsin-EDTA (TLS003703, Worthington) for about 10 min at 28 • C. Subsequently, the digestion of tissue pieces was stopped using 10% fetal bovine serum (FBS) (Gibco) in PBS, and then the cells were filtered through the 40 µm cell strainer followed by 10 min centrifugation at 300 × g, and washed 2-3 times with 1× PBS. Then, 3 ml single cell suspensions was mixed with 5 ml 0.4% dye trypan blue solution to measure the total cell numbers as well as the ratio of live cells using the automated cell counter (Logos Biosystems). Then cells were diluted to the final concentration using 1× PBS with 200 µg/mL bovine albumin (BSA) (NEB, cat # B9000S) for subsequent Drop-seq.
Single Cell Capture, cDNA Library Preparation and Sequencing cDNA Libraries were prepared following the McCarroll Lab Dropseq protocol as previously described (Macosko et al., 2015). Cells were resuspended in PBS with 0.01% BSA and loaded at a concentration of either 150 or 300 cells/µL. Cells were then captured by microfluidic technology in which a single microsphere with a unique molecular identifier (UMI) and a single cell was wrapped in a droplet. The cell capture rate conformed to Poisson distribution, and the best capture rate was about 10%. The cell cleaved in the droplet and released mRNA, which could be captured by microspheres in the same droplet owing to Oligo beads (ChemGenes) that contained the original Drop-seq polyT primer with a VN anchor at the 3 end (TTTTTTTAAGCAGTGGTATCAACGCAGAGTACJJJJJJJJJ JJJNNNNNNNNVTTTTTTTTTTTTTTTTTTTTTTTTTTTTT TVN). The reverse transcription reaction was performed with 25 nt oligo (dT) primer anchored with a 16 nt cell-specific barcode and 10 nt UMIs. Subsequently, the second-strand cDNAs were amplified by PCR. The PCR condition was as follows: 95 • C for 30 s; 11 cycles at 95 • C for 10 s, 55 • C for 30 s, 72 • C for 30 s; followed by 72 • C for 5 min. The RNA-Seq libraries were sequenced on Illumina NextSeq 500.

Classification of Cell Types and Sequencing Analysis
Raw RNA-Sequencing data was analyzed using R package Seurat, v 2.2.0 (Butler et al., 2018). Cell Ranger was performed to output "filtered gene-barcode" count matrix, which was used for multiple downstream analysis. To start, several parameters were used to filter low quality cells according to the expression of mitochondrial genes and cell cycle genes: the total number of expressed genes/cell ranged from 200 to 2500; the total number of UMIs/cell ranged from 300 to 15000; and the percentage of UMIs mapping to mitochondrial genes to total genes was lower than 5%. Moreover, cells were normalized using the "Function NormalizeData" procedure in Seurat. t-SNE was performed through the run_t-SNE function (dims.use = 1:10, max_iter = 2,000) to calculate tSNE plot. Louvain algorithm was performed using Seurat-based workflow. The multiple t-test was performed to obtain the statistical significance of DEGs, and only the genes with significant p-values less than 0.05 with a fold change of log2 transformed FPKM larger than 1.5 were considered as the DEGs. In addition, Function FindMarkers from R package was applied to analyze marker genes between the cells of a cluster and the rest of the cells in the dataset.
To annotate the DEGs function, GO (Gene Ontology) enrichment analysis was implemented by the top GO R packages. Moreover, Gene Set Enrichment Analysis (GSEA 1 ) and were used to perform the KEGG pathway enrichment analysis (Subramanian et al., 2005). The core genes involved in the significantly enriched pathways were performed to generate a heatmap using function heatmap.2 in R package gplots v3.0.1. 1 http://www.webgestalt.org/option.php

In situ Hybridization
In order to verify the scRNA-seq data, three germ cell genes and two somatic genes were validated by in situ hybridization on sections (SISH) as previously reported by Li W. et al. (2017). The ovarian tissues were fixed with 4% paraformaldehade in PBS at 4 • C overnight. After washing with PBS (pH 7.0) three times, the samples were immersed in 30% saccharose-PBS buffer overnight at 4 • C, and then embedded in O.C.T. (Optimal Cutting Temperature, Germany) and sectioned at 10 µm with frozen microtome (Leica, Germany). The cryostat sections were 37 • C for 1 h and then stored at -70 • C. For antisense probe synthesis, T7 RNA polymerase promoter was added to the 5 end of reverse primers using DIG RNA labeling kit (Roche, Germany). Primers sequences for probe synthesis were designed and synthesized according to mRNA sequences in NCBI database (dnd1: XM_018694390.1; nanos3: XM_018664710.1; zar1: XM_018669829.1; hsd17b1: XM_018702509.1; dnajb1: XM_018679787.1) using software Primer 5.0. The primer sequences of these selected genes have given in Supplementary Table 1. Then the RNA probes were synthesized using RNasefree TURBO DNase and purified with LiCl according to the manual of the mMESSAGE mMACHINE kit (cat# AM1340; Ambion). Signal was stained using BCIP/NBT substrates and nuclear was counterstained with propidium iodide (PI). Images were obtained by Zeiss Axio Observer A1 inverted microscope (Leica, Germany).

Conservative and Species Specific Analysis Across Fly, Human, and Asian Sea Bass
Fly and human homologous gene data retrieved from GEO home in National Center for Biotechnology Information Search database (Mahadevaraju et al., 2019, GEO: GSE125948;Zhang et al., 2018, GEO: GSE107746), as well as our data in this study (SRR13619334) were used for the comparative analysis. Housekeeping genes were retrieved as the normalization factor of genes expression from Human Protein Atlas. The common expressed genes among three species were screened using Power Query function in excel and housekeeper genes were picked out from the common genes.

Global Transcriptome Profile of Ovarian Cells in Asian Sea Bass
To construct comprehensive single-cell atlases, single cells were freshly isolated from Asian seabass ovary and the cells were collected using two-step procedure of enzymatic digestion and physical filtering (Valli et al., 2014), then were used for sequencing and for global transcriptome analysis using the Illumina NextSeq 500 platform ( Figure 1A). The morphology of different cells in adult ovarian suspension was shown ( Figure 1B). The large oocytes were separated and deposited at the bottom of the suspension, and the remaining complex contained a number of small oocytes and somatic cells attached to the tissue membrane. Multiple types of cells were mixed together for transcriptome analysis. After standard quality control dataset filters, a total of 1,204 dynamically expressed genes were isolated from 11,094 cells and were retained for downstream analysis (Supplementary Table 2). Saturation analysis revealed that the sequencing depth detection of gene expression was sufficient for further data analysis ( Figure 1C).

Single Cell Transcriptome Profiling of Asian Sea Bass Ovary
To identify distinct cell populations captured and highly variable genes in Asian sea bass ovary, Louvain algorithm was performed using Seurat-based workflow, and the results were visualized via the non-linear dimensionality reduction algorithm t-distributed stochastic neighbor embedding (tSNE) (van der Maaten and Hinton, 2008). As shown in Figure 2A, five clusters of cells were obtained. To display the distinction and similarities across these cell types, the data was filtered using quality control parameters in R package Seurat and the Seurat Violin plots were shown in Supplementary Figure 1. Moreover, the selected genes could clearly be classified the cells into five classes as shown in Figure 2B.
Furthermore, to investigate the biological functions and key pathways of DEGs, Gene Ontology (GO) and the gene set enrichment analysis (GSEA) in Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using WEB-based GEne SeT AnaLysis Toolkit. For GO analysis, "cellular nitrogen compound biosynthetic process, " "cytosol, " and "structural molecule activity" were most enriched GO terms in the category "biological process, " "cellular component, " and "molecular function, " respectively, supporting high levels of cellular activity (Supplementary Figures 2, 3). Importantly, for all five clusters of cells, a total of eight pathways were found, such as biosynthesis of unsaturated fatty acids, ferroptosis, PPAR signaling pathway, progesterone-mediated oocyte maturation, and so on using KEGG analysis ( Figure 3A and Supplementary  Table 4). Notably, the identified functional pathways in this study were reported to be significantly involved in oocytes including the pathway of progesterone-mediated oocyte maturation in described in the previous report (Lutz et al., 2000). In humans, the progesterone-mediated oocyte maturation was significantly overrepresented both in oocytes and granulosa cells of primary stage, suggesting the mediating function in transition from primordial to primary (Lutz et al., 2000). Interestingly, the upregulated genes involved in the eight KEGG pathways including asc15, kras, sqstm1, acox1, and cdk1 were highly expressed in germ cell cluster or granulosa cell cluster (Figure 3B), which indicated that these genes might play significant roles in oocyte maturation or concert germlinegranulosa cells interactions.

Validation of Single Cell Transcriptome Profiling of Asian Sea Bass Ovary
The dendrogram result that cluster 3 was in lineage A and the remaining four cell types were grouped into lineage B suggests that cluster 3 might have distinct gene expression profiles from those of the other four cell types ( Figure 4A). Subsequently, several known cell-type-specific marker genes were used to infer the predominant cell type within each cluster. The expression ddx4, dazl, dnd1, piwil, mre11, rad50 was always considered as a molecular signature for germ cells, and these genes play critical roles in the transition from the primordial to the primary stage during germ cell development or differentiation (Yang et al., 2012;Lotan et al., 2014;Jean et al., 2015;Zhang et al., 2018). Clusters most enriched in kdsr, wt1, and hspg2 were inferred to be granulosa cells as described in previous reports (Zhang et al., 2018). The expression heatmap of these selected genes revealed that most germ cell markers were grouped into cluster 3 while some somatic cell markers were grouped into other four clusters, which confirmed the validity of the tSNE classification ( Figure 4B and Supplementary Table 5). Thus, we defined cluster 3 as germ cells and cluster 5 as granulosa cells, respectively. The other three cell types were unknown cells due to the lack of reported marker genes; thus, we temporarily defined cluster 1, 2, and 4 as ovarian somatic cells.
Moreover, the expression profiles of the known germlinespecific markers including dmc1, dnd1, nanos3, zar1, piwil2, rad50, msh2, sycp3, ddx4, gdf9, figla, and dazl were all     Figure 5A). To further verify the cellular distribution of these marker genes in five clusters, three germ cell related genes, and two somatic cell genes were validated though the in situ hybridization (CISH) on ovarian cryo-sections. The results showed that dnd1 transcripts were detected only in germ cells and absent in surrounding somatic cells, indicating that dnd1 was a germ cell specific gene in Asian sea bass ( Figure 5B). Nanos3 was localized predominantly in the cytoplasm and weakly in the nucleus of oogonia, which was consistent with previous study on zebrafish (Figure 5B; Yamaji et al., 2010). The zar1 transcripts were richly distributed in cytoplasm of oogonia and rarely in somatic cells ( Figure 3B). Cluster 5 was highly expressed marker genes hsd17b1 (Figure 5C), which has been documented to be enriched in granulosa cells in previous reports (Liu et al., 2009). Likewise, in situ hybridization showed that the hsd17b1 transcript were localized to perinuclear cytoplasm of somatic cells inter/peri-oocytes ( Figure 5D). Dnajb1 mRNA signal was mainly detected in the somatic cells surrounding the oocytes and undetectable in oocytes ( Figure 5D).
Subsequently, to further identify and characterize the germ cells and the granulosa cells, the zar1 transcripts enriched cells (zar1 + ) as germ cells and dnajb1 cells (dnajb1 + ) as granulosa cells were sorted out respectively. Totally, 345 cells were sorted out, including 103 cells enriched with dnajb1 expression and 242 cells highly expressing zar1 cells (Supplementary Table 6). Principal component analysis (PCA) performed in R software showed that these cells were obviously clustered into two groups ( Figure 6A). Subsequently the transcriptome profiles of these two cell populations were comparatively analyzed, and it was showed that 44 genes were found to be biased in granulosa (dnajb1 + ) cells, and 596 genes were exclusively expressed in germ cells (zar1 + ), while 124 genes were common in both populations of cells ( Figure 6B and Supplementary Table 6). Especially of the 44 granulosa cell biased genes, 29 genes were new and the first time to be identified in this study (Supplementary Table 7). These differentially expressed genes would be used as biomarkers to trace the interactions between granulosa cells and germ cells during fish ovary developing, and even sex reversing. More importantly, the granulosa cells had also been documented to possess the potency of producing iPSCs, and even oocytes via chemical reprogramming (Mao et al., 2014;Tian et al., 2019). Therefore, the specific cell transcriptome profiles of germ cells and granulosa cells would provide the basic knowledge for further studying the endocrine and fertility functions of the ovarian cells, including granulosa cells and germ cells in fish.

Comparative Analysis of the Ovary Transcriptome Profile Across Phyla Examined
To further investigate the interspecific complexity across fruit fly, human, and Asian sea bass, we compared the ovarian transcriptome profiles of Asian sea bass with previously reported RNA-seq datasets of fruit fly (GEO) home in National Center for Biotechnology Information Search database 2 and human (Zhang et al., 2018). The single-cell data from stages of human oocyte human were merged in order to compare with scRNA-seq data of Asian sea bass ovary and fly oocytes. At last, a total of 26,364 genes expressed in human, 1,709 genes expressed in fly oocytes, and 11,094 genes expressed in Asian sea bass ovarian cells were analyzed. As shown in Figure 7A, there were 151 genes overlapped across fly, human, and Asian sea bass while 5,673 genes shared by human and Asian sea bass, and only 152 genes shared by fly and Asian sea bass. The GO analysis of co-expressed genes was similar with previous investigation (Zhang et al., 2018) that none of these genes were specific or related to folliculogenesis, indicating unique molecular mechanisms underlying the modulation of oocyte development occurred in Asian sea bass ( Figure 7B).
In addition, the conserved and divergent in expression patterns of several oocyte development marker genes such as bmp15, dnd1, patl2, slbp, dmc1, nanos3, piwil1, piwil2, sycp3, ddx4, gdf9, dazl, ybx2 (in mouse, biased expression in adult testis), zar1, and rac1 were defined among human, fish, and fly. Of these genes, only rac1 was detected common in fly, human, and Asian sea bass, the other genes existed only in human and Asian sea bass (Figure 7C), implying that the expression profiles of germ cell genes in fish are more similar with those in human than in fly. Moreover, the common genes shared by Asian sea bass and human also exhibited distinct dynamics expression during oogenesis. As shown in Figure 7C, most genes were relatively expressed at a higher level in human than in Asian sea bass except for nanos3 and zar1. Nanos3 was expressed higher in Asian sea bass than that of key stages of human follicular development in vivo, from primordial to preovulatory stage, while the expression level of zar1 in Asian sea bass was higher than that of human primordial, primary, and preovulatory stages but lower than that of secondary and antral stages. Thus, these findings well documented the conservation and divergence of the shared genes in germ cells across invertebrates and vertebrates 2 https://www.ncbi.nlm.nih.gov/geo and provided a new way to accurately define germ cell and elucidate the mechanism of gonadal development.

DISCUSSION
In the last decade, single-cell transcriptome sequencing has been extensively performed to reveal the molecular signature of generative cells across gametogenesis in invertebrates such as Drosophila (Karaiskos et al., 2017; and mammals such as human, monkey, and mouse (Guo et al., 2013;Knouse et al., 2014;Tsioris et al., 2015;Wang et al., 2020); however, little was reported in lower vertebrates, such as fish, especially the sex reverse species. Both the sex reverse (male-to-female) and ovary development are a complex process in animals, including the Asian sea bass, and thus a full understanding of its regulation requires the integration of multiple data collected from various ovarian cell types. Many studies have shown various interactions or bidirectional communication between female germ cells and surrounding somatic cells (Kloc, 2017), but these kind of data on cellular communications between the somatic cells and oocyte majorly comes from the investigations of mammalian oocytes. Therefore, in this study, we conducted a scRNA-seq analysis to define the signature transcriptome profiles of ovarian cells in Asian seabass adult female, and we identified five cell types, including one type of germ cells and four types of somatic cells. Generally, the transcripts expression profile of ovarian cells in Asian seabass were similar to those discovered in other species. Our study showed that the number of genes or UMIs were more in germ cells than those in somatic cells, which is consistent with the results from fly (Slaidina et al., 2020) and human (Li L. et al., 2017). These findings proved that more transcripts existed in germ cells in gonads. Moreover, in many multicellular organisms, including fish, lots of mRNAs are transcribed during oogenesis and stored in the maternal cytoplasm, providing essential factors for oocyte growth and cell viability, and the mRNAs also are actively translated to regulate embryonic polarity and cell fate during embryogenesis (Olszańska and Borgul, 1993;Seydoux, 1996). It is therefore, critical for precise regulation of the accumulation of proteins implying that a large number of ribosomes must accumulate in the germ cell to translate these messages (Wallace and Selman, 1990). Likewise, the amplification of ribosomal genes allows synthesis of sufficient rRNA to supply the growing oocyte. Besides, most germ cell markers are RNA binding proteins, such as Vasa, Dazl, and Ziwi1 (Kehkooi et al., 2009;Leu and Draper, 2010;Aduma et al., 2019). RNA regulation is, therefore, crucial in the germ line of virtually all organisms and plays a variety of prominent roles. These would be the explanations for the findings of this study: there are many ribosome related genes (rps) were identified in ovarian cells (Supplementary Table 2). Accordingly, the top enriched Go terms include "peptide metabolic process, " "peptide biosynthetic process, " "cellular translation, " "ribonucleoprotein complex, " "RNA binding, " and so on. Furthermore, it is well known that granulosa cells are the most important somatic cells in the ovary, which play important roles in female germ cell development and in the formation of primordial follicles (Kissel et al., 2000). Therefore, it is pivotal to distinguish the granulosa cells and germ cells for further investigations on oogenesis. In our study, we found the zar1 RNA is exclusively expressed in germ cells (Figure 5B), while dnajb1 is specifically expressed in granulosa cells surrounding oocytes ( Figure 5C). Likewise, zygote arrest-1 (zar1) was identified in mice and is considered one of the earliest identified oocytespecific maternal effect gene that functions during MZT (Wu et al., 2003a). Zar1 was also proved to be required for oocyte meiotic maturation by regulating the maternal transcriptome and mRNA translational activation (Rong et al., 2019). Additionally, zar1 is documented to be evolutionarily conserved and expressed in vertebrate ovaries (Wu et al., 2003b). Therefore, we confirmed that the cells expressed zar1 should be germ cells and could be picked out for further analysis. On the other hand, it was reported that the DnaJ-1/Hsp70 chaperone complex regulates border cells (BCs) migration by modulating PVR function during Drosophila oogenesis (Cobreros et al., 2008). Likewise, the germ cells genes were barely detected in the dnajb1 mRNA positive cells surrounding the oocytes, so the dnajb1 + cells were selected as granulosa cells, then the differentially and common genes or pathways were identified between germ cells and granulosa cells in Asian seabass ovary. This would facilitate the investigations on the cellular communications or interactions between somatic cells and germ cells during sex-reversing and/or oogenesis in fish.
Intriguingly, in this study, we also found many cell specific genes, such as hes1 highest in cell cluster C1 in this study, was reported to be highly expressed in granulosa cells and progenitor theca cells in monkey ovary (Zhao et al., 2020). Gdf 9 (highest in cell cluster C3 here) was documented to exhibit high expression level in the oogenesis phase FGCs, which could promote the production of Desert hedgehog (DHH) and Indian hedgehog (IHH) in granulosa cells in monkey ovary (Zhao et al., 2020). The hsd17b1 was found to be enriched in cell cluster C1 and exclusively expressed in most of ovarian somatic cells in our study. The HSD17B1 has been shown to be expressed in fetal human ovaries (Vaskivuo et al., 2005), and the over-expression of human HSD17B1 in mice results in masculinization of the external and internal genitalia of females (Saloniemi et al., 2009); thus, the findings imply that HSD17B1 is likely to contribute in the maintenance of normal sex steroid balance during gonad development. This indicates that hsd17b1 maybe is involved in sex reverse and gonad development in Asian seabass. However, to address this issue, it needs more extensive investigations conducted in the future.
Except accurately distinguishing germ cells from somatic cells, the identifications of the stage-specific genes may provide a valuable clue for the functional studies on gonadal cells or genes, thus illustrating the mechanism behind the gonadal development. The identification of stage-specific marker genes would benefit the understanding of the characterization of cell heterogeneity. In the human ovary, the identification of specific gene signatures of oocytes and granulosa cells in a particular stage may reflect developmental competency and ovarian reserve (Zhang et al., 2018). In adult human testes, the specific marker genes used to identify four undifferentiated spermatogonial (SPG) clusters provide a blueprint of the developing human male germline and supporting somatic cells (Sohni et al., 2019). In the mouse, dynamic changes in the transcriptomes of spermatogenic cells were observed based on the expression level of key cell type-specific markers (Grive et al., 2019;Suzuki et al., 2019). So far, numerous efforts made to identify biomarkers are mainly focused on mammals; few marker genes have been identified in lower vertebrate such as fish, and their specificity is still largely unclear. In this study, we identified germ cell and somatic cell populations based on several reported candidate biomarkers (Figure 2). We also detected a set of new candidate genes (Supplementary Table 2), which would provide a basic and valuable index for demonstrating the ovarian development, and also are critical for developing related bio-techniques, such as isolating, culturing, and identifying the germ cells in fish or other vertebrates. We uncovered several functional signaling pathways that significantly represented those associated with oocyte meiosis, including progesterone-mediated oocyte maturation (Figure 3). Genes involved in these functional pathways were highly expressed in germ cells (cluster 3) or granulosa cells (cluster 5), indicating that these functional pathways may concert germline-somatic interactions. Hsp90ab1 gene in progesterone-mediated oocyte maturation pathway has been demonstrated to be required for male fertility, sex determination, and metabolic homeostasis through interaction with Kdm3a (Ioannis et al., 2014). Moreover, ccna2, ccnb1, ccnb3, cdk1, cdc26, kras, and plk1 were involved in multiple signal pathways. The connection and interaction of genes indicated that the ovarian development was a complex process mediated by multiple factors and pathways.
Understanding the transition from mitosis to meiosis in early oogenesis is essential for investigating ovary development and the application of induced gametes. In previous reports, it has been found that the oocyte transcriptome is variable across animal kingdom (Biase, 2017;Gahurova et al., 2017). For example, strong transcriptional activity and considerable variability were uncovered in both human and mouse oocytes (Sylvestre et al., 2013). In the present study, we compared gene expression dynamics of invertebrate fruit flies, lower vertebrate Asian sea bass, and human ovarian cell transcriptome, and we found distinction in gene expression profiles across these three species with less difference between Asian sea bass and humans than that between Asian sea bass and fruit flies. On the other hand, a certain degree of similarity in gene expression was also detected among these three species, and the further function analysis of the conserved genes may provide key insights in reproductive biology and fertility in metazoans. Therefore, our research on ovarian cell transcriptome of Asian sea bass can provide basis for further studies on the molecular mechanisms associated with ovarian development in other fish, even in human.
In summary, this study is the first comprehensive single-cell transcriptomic atlas of fish ovary, especially a hermaphrodite fish. The cell identities and cell-type-specific gene signatures in the fish ovary were defined; likewise, the findings lay a foundation for the future investigations on the identification and isolation of ovarian cells in fish. More importantly, the findings offer new insights into understanding the conservation and divergence in the molecular mechanisms underlying ovarian cells development across phyla.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI Sequence Read Archive database under the accession number of SRR13619334 (Biosample SAMN17718760).

ETHICS STATEMENT
The animal study was reviewed and approved by the Laboratory Animal Ethics Committee Pearl River Fisheries Research Institute, CAFS.

AUTHOR CONTRIBUTIONS
HX, XL, and XZ conceived the project and designed the experiments. WL, XL, and YY conducted the experiments. KC, HY, and YL analyzed the results. HX and XL wrote the manuscript. All the authors read and approved the final manuscript.