The Comparative Survey of Coordinated Regulation of Steroidogenic Pathway in Japanese Flounder (Paralichthys olivaceus) and Chinese Tongue Sole (Cynoglossus semilaevis)

Steroidogenesis controls the conversion of cholesterol into steroid hormones through the complex cascade reaction of various enzymes, which play essential roles in sexual differentiation and gonadal development in vertebrates, including teleosts. Japanese flounder (Paralichthys olivaceus) and Chinese tongue sole (Cynoglossus semilaevis) are important marine cultured fishes in China and have remarkable sexual dimorphism with bigger females and sex reversal scenarios from female to neo-male. Several steroidogenic genes have been analyzed individually in the two species, but there is a lack of information on the coordinated interaction of steroidogenic gene regulation. Therefore, in this study, through genomic and transcriptomic analysis, 39 and 42 steroidogenic genes were systematically characterized in P. olivaceus and C. semilaevis genomes, respectively. Phylogenetic and synteny analysis suggested a teleost specific genome duplication origin for cyp19a1a/cyp19a1b, hsd17b12a/hsd17b12b, ara/arb and esr2a/esr2b but not for star/star2 and cyp17a1/cyp17a2. Comparative transcriptome analysis revealed conserved expression patterns for steroidogenic genes in P. olivaceus and C. smilaevis gonads; star/star2, cyp11a/cyp11c, cyp17a1/cyp17a2, cyp21a, hsd3b1, hsd11b and hsd20b were strongly expressed in testis, while cyp19a1a and hsd17b genes were highly expressed in ovaries. Only a few genes were differentially expressed between male and neo-male testis of both P. olivaceus and C. semilaevis, and even fewer genes were differentially regulated in the brains of both species. Network analysis indicated that cyp11c, cyp17a1 and hsd3b1 actively interacted with other steroidogenic genes in P. olivaceus and C. semilaevis, and may play a more sophisticated role in the steroid hormone biosynthesis cascade. The coordinated interaction of steroidogenic genes provided comprehensive insights into steroidogenic pathway regulation with a global biological impact, as well as sexual development in teleost species.


Introduction
Steroidogenesis is the conversion of cholesterol into biologically active steroid hormones by the complex cascade reaction of various enzymes, which plays essential roles during many physiological processes including sexual differentiation, gonadal development, growth and maturation [1]. Steroidogenesis occurs primarily in tissues, such as the gonads, interrenal gland and brain, which is, in general, controlled by the hypothalamus-pituitary-interrenal

Genomic Landscape, Functional Domain and Phylogeny of Steroidogenic Genes in P. olivaceus and C. semilaevis
Through a genome-wide screen, 39 and 42 genes participating in steroid biosynthesis were identified in the P. olivaceus and C. semilaevis genomes, respectively ( Figure 1 and Table S1), including star, cyp, hsd genes, as well as steroid receptors genes esr and ar ( Figure 1 and Table S2). The copy number of these steroidogenic genes between tetrapod and teleost lineages was generally conserved with several gene duplication events observed in the teleost lineage (Table S1), most likely originating from the teleost specific genome duplication (TSGD).

Genomic Landscape, Functional Domain and Phylogeny of Steroidogenic Genes in P. olivaceus and C. semilaevis
Through a genome-wide screen, 39 and 42 genes participating in steroid biosynthesis were identified in the P. olivaceus and C. semilaevis genomes, respectively ( Figure 1 and Table S1), including star, cyp, hsd genes, as well as steroid receptors genes esr and ar ( Figure  1 and Table S2). The copy number of these steroidogenic genes between tetrapod and teleost lineages was generally conserved with several gene duplication events observed in the teleost lineage (Table S1), most likely originating from the teleost specific genome duplication (TSGD). Specifically, 16 and 17 cyp genes were obtained in the P. olivaceus and C. semilaevis genomes, respectively. There were comparable numbers of cyp1, cyp11, cyp21, cyp26 and cyp27 genes between tetrapods and teleosts, whereas two cyp17s and two cyp19s were identified in teleosts but only one copy was found in tetrapods, indicating their putative TSGD origin in teleost lineages ( Figure 1 and Table S1). Five conserved motifs (motif [1][2][3][4][5] were shared among P. olivaceus cyp genes with motif 6 only present in cyp1s, cyp17s and cyp21, which was consistent with their phylogenetic relationships (Figure 2a, b). Moreover, the vertebrate cyp genes were clustered into seven clades, corresponding to the cyp1, cyp11, cyp17, cyp19, cyp21, cyp26 and cyp27 families (Figure 2b). The cyp17 clade was further divided into two branches with the teleost cyp17a1s clustered with orthologous of tetrapods, while teleost cyp17a2s were further clustered together (Figure 2b). The cyp19a1 clade was also grouped into two branches, which corresponded to cyp19a1 of tetrapods and cyp19a1a/cyp19a1b of teleosts, indicating their TSGD origin in teleost lineages. Specifically, 16 and 17 cyp genes were obtained in the P. olivaceus and C. semilaevis genomes, respectively. There were comparable numbers of cyp1, cyp11, cyp21, cyp26 and cyp27 genes between tetrapods and teleosts, whereas two cyp17s and two cyp19s were identified in teleosts but only one copy was found in tetrapods, indicating their putative TSGD origin in teleost lineages ( Figure 1 and Table S1). Five conserved motifs (motif 1-5) were shared among P. olivaceus cyp genes with motif 6 only present in cyp1s, cyp17s and cyp21, which was consistent with their phylogenetic relationships (Figure 2a,b). Moreover, the vertebrate cyp genes were clustered into seven clades, corresponding to the cyp1, cyp11, cyp17, cyp19, cyp21, cyp26 and cyp27 families (Figure 2b). The cyp17 clade was further divided into two branches with the teleost cyp17a1s clustered with orthologous of tetrapods, while teleost cyp17a2s were further clustered together (Figure 2b). The cyp19a1 clade was also grouped into two branches, which corresponded to cyp19a1 of tetrapods and cyp19a1a/cyp19a1b of teleosts, indicating their TSGD origin in teleost lineages. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 4 of 17 Figure 2. Characterization of steroidogenic proteins. The conserved motif and structure of (a) cyp, (c) hsd, (e) star, esr and ar in P. olivaceus, and the phylogeny of (b) cyp, (d) hsd, (f) star, esr and ar in selected vertebrates. The phylogenic trees were constructed with a bootstrap of 1,000 times. The subfamilies of cyp and hsd were labeled with different background colors, and the P. olivaceus proteins were highlighted in red.

Figure 2.
Characterization of steroidogenic proteins. The conserved motif and structure of (a) cyp, (c) hsd, (e) star, esr and ar in P. olivaceus, and the phylogeny of (b) cyp, (d) hsd, (f) star, esr and ar in selected vertebrates. The phylogenic trees were constructed with a bootstrap of 1000 times. The sub-families of cyp and hsd were labeled with different background colors, and the P. olivaceus proteins were highlighted in red.
Hsd3b, hsd11b, hsd17b and hsd20b were pivotal members of the hsd genes. In total, 16 hsd genes were obtained in the P. olivaceus and C. semilaevis genomes, respectively, which were generally conserved between tetrapods and teleosts, with hsd20b only present in teleosts but not in tetrapods ( Figure 1). Two or three hsd3b genes were obtained in most selected vertebrates as the conserved hsd3b1 and hsd3b7, while hsd3b1, -2, -3, -4, -5, -6 only existed in mice (Table S2). Hsd17b was the largest hsd family, with 11 and 10 hsd17b genes from P. olivaceus and C. semilaevis, respectively ( Figure 1). The functional domain of P. olivaceus hsd genes was less conserved with 12 motifs detected ( Figure 2c). Moreover, hsd17b12a and hsd17b12b were clustered into two branches in teleosts and further clustered with tetrapod hsd17b12s, suggesting a genome duplication origin in teleosts from the ancestor hsd17b12 ( Figure 2d).
In addition, a single copy of star was present in tetrapods while two copies were identified in most teleost species as star and star2 ( Figure 1). Interestingly, three copies of star existed in C. semilaevis (Table S2), with star-2 and star-3 located on Z and W chromosomes of C. semilaevis, respectively. With regard to esr, esr1 was conserved as a single copy among selected vertebrates, whereas there were two esr2s, esr2a and esr2b, in most teleost species (Table S2). This was also found in ar genes, with ara and arb present in the teleost lineages (Table S2). In P. olivaceus, the functional motifs of esr, ar and star genes were generally conserved. As sex hormone receptors, esr and ar shared seven motifs, with only motif 5 and motif 7 being specific to their C-terminals, respectively ( Figure 2e). Moreover, star and star2 were split into two clades, which was similar to the findings of Yu et al. [32]. Esrs were also clustered into two clades, as esr1 and esr2, with esr2 further divided into two branches (Figure 2f), as esr2 in tetrapods and esr2a and esr2b in teleosts.

The Duplication of Steroidogenic Genes in the Teleost Lineages
Teleost fishes and humans show remarkable conservation in many developmental and physiological aspects, including the endocrine system and steroid hormone related processes [6]. Here, to address the impact of genome duplication in teleosts on steroid hormone biosynthesis, seven pairs of steroidogenic genes, including star/star2, cyp11a/cyp11c, cyp17a1/cyp17a2, cyp19a1a/cyp19a1b, hsd17b12a/hsd17b12b, esr2a/esr2b and ara/arb, were employed for the synteny analysis.
Firstly, there was conserved synteny between cyp19a1a and cyp19a1b in vertebrates. For example, conservative upstream (gldn, dmxl2) and downstream (tnfaip8l3) genes were found in both cyp19a1 and cyp19a1a, while ap4e1 was adjacent to cyp19a1 and cyp19a1b ( Figure 3a). It was noteworthy that cyp19a1b was absent in fugu, but the conserved neighboring genes and chromosomes can still be identified ( Figure 3a). This indicated that cyp19a1a and cyp19a1b originated from TSGD, which supported the previous conclusion [33]. In addition, both hsd17b12a and hsd17b12b in teleosts were found to share conserved adjacent genes with hsd17b12, such as alkbh3 and api5 between hsd17b12 of higher vertebrates and hsd17b12a of teleosts, and ttc17 and lrrc4c between hsd17b12 and hsd17b12b ( Figure 3b). Consistent with the phylogeny, esr2 in tetrapods shared the common gene syne2 with teleost esr2a and esr2b (Figure 3c). Moreover, esr2 and esr2b were more conservative with neighboring genes wdr89, ppp2r5e and gphb5 (Figure 3c). Similarly, ar and ara contained the conserved neighboring gene efnb1, while ophn1 was the common gene between ar and arb (Figure 3d), consistent with their phylogeny (Figure 2f). Together with the results of phylogeny and synteny, it suggested that cyp19a1a/1b, hsd17b12a/b, ara/b and esr2a/2b in teleosts indeed originated from the TSGD. On the other hand, the synteny of star between tetrapods and teleosts was generally conserved but not with star2 of teleosts ( Figure S1a). For example, grk5, lsm1 and bag4 were conserved downstream genes of star, with grk5 absent in humans, mice and chickens. Upstream genes of star were less conserved, with ash2l existing only in tetrapods and spotted gar ( Figure S1a). This indicated that although two star genes were present in teleosts, star2 probably originated independently in the teleost lineage during evolution. This also can be seen in cyp11s and cyp17as ( Figure S1b, c). For instance, no conserved upstream and downstream genes were identified between cyp17a1 and cyp17a2 in teleosts, proving that they did not originate from the TSGD (Figure S1c), although cyp17a1 had conserved downstream genes between tetrapods and teleosts, such as borcs7 and nt5c2. Moreover, the phylogeny of cyp17a1 between tetrapods and teleosts was closer than that of cyp17a1 and On the other hand, the synteny of star between tetrapods and teleosts was generally conserved but not with star2 of teleosts ( Figure S1a). For example, grk5, lsm1 and bag4 were conserved downstream genes of star, with grk5 absent in humans, mice and chickens. Upstream genes of star were less conserved, with ash2l existing only in tetrapods and spotted gar ( Figure S1a). This indicated that although two star genes were present in teleosts, star2 probably originated independently in the teleost lineage during evolution. This also can be seen in cyp11s and cyp17as ( Figure S1b,c). For instance, no conserved upstream and downstream genes were identified between cyp17a1 and cyp17a2 in teleosts, proving that they did not originate from the TSGD (Figure S1c), although cyp17a1 had conserved downstream genes between tetrapods and teleosts, such as borcs7 and nt5c2. Moreover, the phylogeny of cyp17a1 between tetrapods and teleosts was closer than that of cyp17a1 and cyp17a2 (Figure 2b), indicating that cyp17a1 and cyp17a2 in teleost species were probably gained independently during evolution. Therefore, it revealed that star/star2, cyp11a/11c and cyp17a1/17a2 did not originate from TSGD as previously assumed [6].

Evolutionary Dynamics of Steroidogenic Genes in Teleost Lineages
To further understand the evolutionary dynamics of steroidogenic pathways in teleost lineages, the ratio of non-synonymous substitution to synonymous substitution (ω, dN/dS) was analyzed to detect the possibility of selection. Site model (SM) tests showed generally conserved ω values (0-0.4) among teleost steroidogenic genes (Table S3 and Figure S2a), indicating their functional conservation during evolution. Specifically, hsd20b2 represented the highest ω values (0.3941), while cyp26b1 showed the lowest (0.0291) ( Figure S2a). There were also diverse ω values in closely related genes, such as between star and star2, cyp26a1/b1/c1, hsd11b1L and hsd11b2, and hsd17b genes ( Figure S2a). In branch model (BM) tests, the ω values of nine genes were significantly different between P. olivaceus and other teleost species ( Figure S2b), among which seven genes of cyp11c, cyp21, cyp27c1, hsd11b1L, hsd17b4, hsd17b8 and hsd17b14 exhibited higher ω values in P. olivaceus than in other teleost species, indicating greater selection pressure in P. olivaceus ( Figure S2b). Moreover, among the putatively duplicated steroidogenic genes, BM tests indicated varied evolutionary patterns. For instance, in teleost lineages, star2, cyp17a2 and cyp19a1b showed a faster evolutionary rate (ω1) than their paralogues (ω0) (Figure 4), while other TSGD originated genes (hsd17b12a/b, ara/b and esr2a/2b) do not represent elevated evolution ( Figure 4), indicating possible divergence independent from the TSGD.

Evolutionary Dynamics of Steroidogenic Genes in Teleost Lineages
To further understand the evolutionary dynamics of steroidogenic pathways in teleost lineages, the ratio of non-synonymous substitution to synonymous substitution (ω, dN/dS) was analyzed to detect the possibility of selection. Site model (SM) tests showed generally conserved ω values (0-0.4) among teleost steroidogenic genes (Table S3 and  ). There were also diverse ω values in closely related genes, such as between star and star2, cyp26a1/b1/c1, hsd11b1L and hsd11b2, and hsd17b genes ( Figure S2a). In branch model (BM) tests, the ω values of nine genes were significantly different between P. olivaceus and other teleost species ( Figure S2b), among which seven genes of cyp11c, cyp21, cyp27c1, hsd11b1L, hsd17b4, hsd17b8 and hsd17b14 exhibited higher ω values in P. olivaceus than in other teleost species, indicating greater selection pressure in P. olivaceus ( Figure  S2b). Moreover, among the putatively duplicated steroidogenic genes, BM tests indicated varied evolutionary patterns. For instance, in teleost lineages, star2, cyp17a2 and cyp19a1b showed a faster evolutionary rate (ω1) than their paralogues (ω0) (Figure 4), while other TSGD originated genes (hsd17b12a/b, ara/b and esr2a/2b) do not represent elevated evolution ( Figure 4), indicating possible divergence independent from the TSGD.

Expression Profiles of Steroidogenic Genes in P. olivaceus and C. semilaevis
The steroidogenesis during sexual development in teleosts is essentially regulated by differential expression of several steroidogenic enzymes. Through comparative transcriptome analysis of brains and gonads, most steroidogenic genes represented conserved expression patterns between P. olivaceus and C. semilaevis (Table S4). In gonad tissues, the most differentially expressed genes (DEGs) were obtained between the ovaries and testis (male/neo-male), while only a few DEGs were found between the normal male and sex reversed neo-male testes (Figures 5 and S3). For example, in the two gynogenetic P. olivaceus gonad transcriptomes, cyp11a, cyp11c, cyp17a1 and hsd3b1 had testis biased expression (Figure 5a, b), while cyp19a1a, cyp26a1, hsd17b1 and hsd17b12a had ovary biased expression (Figure 5a, b). In addition, the expression of steroidogenic genes between the gonads of the wide type and gynogenetic P. olivaceus was also similar (Figure 5a), with only cyp26a1 being highly expressed in the gynogenetic female ovaries but not in the wild type P. olivaceus gonads (Figure 5a), indicating consistent steroid biosynthesis functions between wild type and gynogenetic individuals. Similarly, in the C. semilaevis gonads, there were

Expression Profiles of Steroidogenic Genes in P. olivaceus and C. semilaevis
The steroidogenesis during sexual development in teleosts is essentially regulated by differential expression of several steroidogenic enzymes. Through comparative transcriptome analysis of brains and gonads, most steroidogenic genes represented conserved expression patterns between P. olivaceus and C. semilaevis (Table S4). In gonad tissues, the most differentially expressed genes (DEGs) were obtained between the ovaries and testis (male/neo-male), while only a few DEGs were found between the normal male and sex reversed neo-male testes (Figures 5 and S3). For example, in the two gynogenetic P. olivaceus gonad transcriptomes, cyp11a, cyp11c, cyp17a1 and hsd3b1 had testis biased expression (Figure 5a,b), while cyp19a1a, cyp26a1, hsd17b1 and hsd17b12a had ovary biased expression (Figure 5a,b). In addition, the expression of steroidogenic genes between the gonads of the wide type and gynogenetic P. olivaceus was also similar (Figure 5a), with only cyp26a1 being highly expressed in the gynogenetic female ovaries but not in the wild type P. olivaceus gonads (Figure 5a), indicating consistent steroid biosynthesis functions between wild type and gynogenetic individuals. Similarly, in the C. semilaevis gonads, there were also conserved expression patterns of steroidogenic genes. For example, star-3, cyp11a/c, cyp17a1/a2, cyp21a, cyp26a1, hsd3b1/b7, hsd11b2, hsd20b2, esr1 and esr2a all had testis biased expression in the two transcriptomes, while only hsd17b12a had ovary biased expression (Figure 5c,d). Interestingly, there were more testis biased than ovary biased steroidogenic genes in both the P. olivaceus and C. semilaevis gonads, which was in line with our previous finding [34], suggesting significant male sex steroid hormone activities in both the P. olivaceus and C. semilaevis testis. Furthermore, there was less than expected DEGs between the normal male and sex reversed neo-male testis of both P. olivaceus and C. semilaevis. For example, only cyp19a1b was differentially expressed between the neo-male and male testis of P. olivaceus, and no DEGs were present between the neo-male and male testis of C. semilaevis (Figure 5a,d). This was probably due to the similar steroid hormone biosynthesis process between the mature male and neo-male testis, which could both develop functional sperms. Therefore, other functional pathways may warrant further investigation to ascertain specific differences in testis development and spermatogenesis between male and neo-male P. olivaceus and C. semilaevis. In addition, only a few steroidogenic genes were differentially expressed between the male and female brains of P. olivaceus and C. semilaevis. For example, four genes were differentially expressed between the brains of male and female gynogenetic P. olivaceus, with cyp21a and hsd3b1 being male brain biased, and cyp19a1a and hsd17b1 being female brain biased ( Figure 5b); only star-3 was differentially expressed between the brains of male and female C. semilaevis (Figure 5c).
To further understand the function of steroidogenic genes, GO and KEGG enrichment of their homologous in model organisms (human and zebrafish) were analyzed with various biological processes (BPs), molecular functions (MFs) and cellular components (CCs) ( Figure S4). Specifically, the steroidogenic genes were enriched in five BPs, including steroid biosynthetic process, estrogen biosynthetic process, response to estrogen, glucocorticoid biosynthetic process and sterol metabolic process ( Figure S4a); five MFs contained estradiol 17-beta-dehydrogenase activity, testosterone dehydrogenase (NAD+) activity, steroid binding, aromatase activity and steroid hydroxylase activity ( Figure S4b), as well as one CC of the endoplasmic reticulum membrane ( Figure S4c). Moreover, the steroidogenic genes were enriched in five KEGG pathways, including steroid hormone biosynthesis, metabolic pathways, retinol metabolism, ovarian steroidogenesis, aldosterone synthesis and secretion ( Figure S4d). Therefore, these genes were supported to participate in various processes of steroid hormone biosynthesis. For example, cyp19a1a and several hsd17b (1,2,8,11,12) genes were mainly associated with the estrogen biosynthetic process and estradiol 17-beta-dehydrogenase activity ( Figure S4a,b), which was in agreement with their ovary biased expression pattern ( Figure 5). Moreover, cyp11a/c, cyp21a, hsd3b1/b2 and hsd11b1/b2 genes participated in both glucocorticoid biosynthetic processes and aldosterone synthesis and secretion ( Figure S4a

Coordinated Interaction of Steroidogenic Pathway in P. olivaceus and C. semilaevis
To investigate the possible cascade regulation of steroidogenic genes, the STRING database was employed to show that the steroid biosynthesis proteins could interact directly with each other in model organisms ( Figure S5). For example, cyp1a, cyp11b/c, cyp17a1, cyp19a1, hsd3b1 and several hsd17b genes were highly interactive with other steroidogenic genes in zebrafish and humans ( Figure S5). In addition, with PCC > 0.8 as highly correlation, interaction networks revealed a high level of coordinated regulation among steroidogenic genes of P. olivaceus and C. semilaevis, with more positive correlation than negative correlation (Table S5). For instance, in P. olivaceus, cyp11c, cyp17a1, cyp19a1a, hsd3b1, hsd20b and hsd17b genes appeared to interact greatly with other steroidogenic genes (Figure 6a,b), while cyp11a/c, cyp17a1/a2, cyp21a, hsd3b1, hsd11b2 and hsd20b2 were the most interactive genes in C. semilaevis (Figure 6c,d). Most of these genes represented testis biased expression ( Figure 5), indicating their possible essential coordinated roles in testicular development and spermatogenesis. Moreover, in P. olivaceus and C. semilaevis, the sex hormone receptor genes, esr and ar, also represented high connectivity with other genes (Figure 6), which was in accordance with their transcription factor function. For example, in the steroidogenic pathway, both esr and ar may have direct interactions with cyp19a1 genes, and therefore, participate in the regulation of the transcriptional activity of cyp19a1 [31].

Coordinated Interaction of Steroidogenic Pathway in P. olivaceus and C. semilaevis
To investigate the possible cascade regulation of steroidogenic genes, the STRING database was employed to show that the steroid biosynthesis proteins could interact directly with each other in model organisms ( Figure S5). For example, cyp1a, cyp11b/c, cyp17a1, cyp19a1, hsd3b1 and several hsd17b genes were highly interactive with other steroidogenic genes in zebrafish and humans ( Figure S5). In addition, with PCC > 0.8 as highly correlation, interaction networks revealed a high level of coordinated regulation among steroidogenic genes of P. olivaceus and C. semilaevis, with more positive correlation than negative correlation (Table S5). For instance, in P. olivaceus, cyp11c, cyp17a1, cyp19a1a, hsd3b1, hsd20b and hsd17b genes appeared to interact greatly with other steroidogenic genes (Figure 6a,b), while cyp11a/c, cyp17a1/a2, cyp21a, hsd3b1, hsd11b2 and hsd20b2 were the most interactive genes in C. semilaevis (Figure 6c, d). Most of these genes represented testis biased expression ( Figure 5), indicating their possible essential coordinated roles in testicular development and spermatogenesis. Moreover, in P. olivaceus and C. semilaevis, the sex hormone receptor genes, esr and ar, also represented high connectivity with other genes (Figure 6), which was in accordance with their transcription factor function. For example, in the steroidogenic pathway, both esr and ar may have direct interactions with cyp19a1 genes, and therefore, participate in the regulation of the transcriptional activity of cyp19a1 [31]. In the sophisticated reaction cascade of steroid hormone biosynthesis in the gonads, star assists the transport of precursor cholesterol from the outer to the inner mitochondrial membrane (Figure 7) [4]. In P. olivaceus and C. semilaevis, star (star-2/3) and star2 were highly expressed in the testis of both species (Figures 5 and 7), indicating their functional conservation in the production of male-specific steroid hormones. In Nile tilapia (Oreochromis niloticus), androgen generation during spermatogenesis may require the involvement of star, which was proven to be primarily expressed in the testis [32]. Interestingly, even with star-1 not expressed in the gonads of C. semilaevis, star-2 and star-3 were highly expressed in the testis of C. semilaevis (Figure 5c, d), suggesting their diverged function in C. semilaevis testis development and spermatogenesis. In mitochondria, cyp11a is responsible for subsequently catalyzing the side-chain hydroxylation and cleavage of cholesterol In the sophisticated reaction cascade of steroid hormone biosynthesis in the gonads, star assists the transport of precursor cholesterol from the outer to the inner mitochondrial membrane (Figure 7) [4]. In P. olivaceus and C. semilaevis, star (star-2/3) and star2 were highly expressed in the testis of both species (Figures 5 and 7), indicating their functional conservation in the production of male-specific steroid hormones. In Nile tilapia (Oreochromis niloticus), androgen generation during spermatogenesis may require the involvement of star, which was proven to be primarily expressed in the testis [32]. Interestingly, even with star-1 not expressed in the gonads of C. semilaevis, star-2 and star-3 were highly expressed in the testis of C. semilaevis (Figure 5c,d), suggesting their diverged function in C. semilaevis testis development and spermatogenesis. In mitochondria, cyp11a is responsible for subsequently catalyzing the side-chain hydroxylation and cleavage of cholesterol to pregnenolone (Figure 7) [35], which together with star for cholesterol mobilization, controls the rate of steroid hormone biosynthesis [22]. Cyp11a was greatly expressed in both the male and neo-male testes of P. olivaceus and C. semilaevis (Figures 5 and 7), and it interacted with star in C. semilaevis network (Figure 6c), indicating that the steroidogenesis is officially initiated under the action of star and cyp11a, which may be more closely related to testicular development in teleost species [36].
to pregnenolone (Figure 7) [35], which together with star for cholesterol mobilization, controls the rate of steroid hormone biosynthesis [22]. Cyp11a was greatly expressed in both the male and neo-male testes of P. olivaceus and C. semilaevis (Figures 5 and 7), and it interacted with star in C. semilaevis network (Figure 6c), indicating that the steroidogenesis is officially initiated under the action of star and cyp11a, which may be more closely related to testicular development in teleost species [36]. Cyp17a and hsd3b are important enzymes as the second step of steroidogenesis involving the conversion of pregnenolone into 17α-OH-pregnenolone and progesterone, and later to androstenedione [12,14,37] (Figure 7). Cyp17a1 and cyp17a2 were reported to be highly expressed in the P. olivaceus testis and ovary [30], and cyp17a1 and hsd3b1/b7 were all strongly expressed in the testis of both P. olivaceus and C. semilaevis ( Figure 5). Here, in the networks, cyp17a1/a2 and hsd3b1 were closely associated with each other (Figure 6), which corresponded to their catalytic function of transforming pregnenolone into androstenedione (Figure 7). Moreover, cyp17a1/a2 and hsd3b1/b7 also strongly interacted with testis biased cyp11a/c, cyp21a, hsd20b, as well as negatively correlated with ovary biased cyp19a1a and hsd17b genes (Figure 6), indicating their essential regulatory interaction in the steroidogenesis cascade. For example, hsd20b has been demonstrated to be associated with 17α, and 20β-DP synthesis [38] as well as testicular development [39]. Here hsd20b was dominant in testicular expression in both P. olivaceus and C. semilaevis ( Figure  5) and has a similar interaction pattern as hsd3b1 with other steroidogenic genes ( Figure  6), corresponding to the sequential interaction of hsd3b1 and hsd20b along with the biosynthesis of the C21 steroid (Figure 7). Cyp17a and hsd3b are important enzymes as the second step of steroidogenesis involving the conversion of pregnenolone into 17α-OH-pregnenolone and progesterone, and later to androstenedione [12,14,37] (Figure 7). Cyp17a1 and cyp17a2 were reported to be highly expressed in the P. olivaceus testis and ovary [30], and cyp17a1 and hsd3b1/b7 were all strongly expressed in the testis of both P. olivaceus and C. semilaevis ( Figure 5). Here, in the networks, cyp17a1/a2 and hsd3b1 were closely associated with each other (Figure 6), which corresponded to their catalytic function of transforming pregnenolone into androstenedione (Figure 7). Moreover, cyp17a1/a2 and hsd3b1/b7 also strongly interacted with testis biased cyp11a/c, cyp21a, hsd20b, as well as negatively correlated with ovary biased cyp19a1a and hsd17b genes (Figure 6), indicating their essential regulatory interaction in the steroidogenesis cascade. For example, hsd20b has been demonstrated to be associated with 17α, and 20β-DP synthesis [38] as well as testicular development [39]. Here hsd20b was dominant in testicular expression in both P. olivaceus and C. semilaevis ( Figure 5) and has a similar interaction pattern as hsd3b1 with other steroidogenic genes (Figure 6), corresponding to the sequential interaction of hsd3b1 and hsd20b along with the biosynthesis of the C21 steroid ( Figure 7).
Subsequently, cyp11c and hsd11b participate in the production of 11-ketotestosterone (11-KT) precursors as well as the biosynthesis of cortisol in teleosts ( Figure 7); therefore, they could function essentially in the crosstalk between the HPG and HPI axes [40,41]. In addition, steroid 21-hydroxylase (cyp21a) has a crucial role in the formation of 11deoxycorticosterone and 11-deoxycortisol [6]. Together with cyp11c, hsd11b was proven to synthesize 11-KT from testosterone (T) and interact with cyp21a to promote the formation of cortisol and corticosterone (Figure 7) [42], which was consistent with the reaction network results in P. olivaceus and C. semilaevis that cyp11c was strongly associated with hsd11b and cyp21a, respectively ( Figure 6). All these three genes were strongly expressed in the male and neo-male testis of P. olivaceus and C. semilaevis ( Figure 5), indicating their essential roles in androgen related hormones biosynthesis. Interestingly, hsd11b1la was highly expressed in the P. olivaceus testis but not in the C. semilaevis testis, while hsd11b2 was greatly expressed in the C. semilaevis testis but not in the P. olivaceus testis (Figure 5), suggesting their putatively lineage-specific function in different teleost species.
On the other side, aromatase (cyp19a1a) catalyzes the production of aromatic C18 estrogens from C19 androgens, for instance, from androstenedione into estrone (E1) and from T into estradiol (E2) (Figure 7) [28,43]. The expression of cyp19a1a was higher in the ovary than in the testis, while cyp19a1b was notably highly expressed in the brain than in the gonads as well as highly expressed in neo-male brains of P. olivaceus and C. semilaevis (Table S4). A similar cyp19a1 expression pattern was also reported in other teleost species, such as in Clupeocephala and zebrafish [44,45], which suggested that cyp19a1b is relevant to brain development and repair in teleost lineages [46]. Together with cyp19a1a, hsd17b1 is responsible for the generation of E2 [47], and the network results were consistent with this collaboration, showing a strong correlation between cyp19a1a and hsd17b genes (Figure 6a,c,d). Hsd17b1 was proven to be highly expressed in the ovary of P. olivaceus [24], but over-expressed in the testis of C. semilaevis ( Figure 5). Hsd17b3 is involved in T production [20], with the expression of hsd17b3 higher in the testis of O. niloticus and P. olivaceus [24,48]. In this study, hsd17b3 was not differentially expressed in either P. olivaceus or C. semilaevis gonads ( Figure 5), which requires further experimental verification. Moreover, hsd17b12a and hsd17b12b can accelerate the conversion of E1 into E2 and promote estrogen synthesis [49]. Only hsd17b12a but not hsd17b12b was highly expressed in the ovary of P. olivaceus and C. semilaevis ( Figure 5), suggesting their functional divergence. This was in line with the previous report that hsd17b12 was distinctly expressed in the P. olivaceus ovary [24]. All in all, these interactive networks indicate the cascading effect of key steroidogenic enzymes and warrant further investigation with a global biological impact in teleost lineages.
Multiple Em for Motif Elicitation (MEME, https://meme-suite.org/meme (accessed on 7 July 2021)) program was applied to evaluate the conserved motifs with the multiple sequence alignment of the aa sequences, with the parameters of any number of repetitions, optimum width of motifs from 6 to 200 with 6 motifs for cyp genes, 12 motifs for hsd genes, 12 motifs for star, ar and esr, respectively. The Gene Structure Display Sever (GSDS 2.0, http://gsds.gao-lab.org (accessed on 7 July 2021)) was used to visualize the exon-intron structure of steroidogenic genes. The gene structure and conserved motifs were visualized using TBtools V1.09 [50].

Phylogenetic Analysis of Steroidogenic Genes
The multiple sequence alignment was performed using MUSCLE [51] with the aa sequences of cyp, hsd, star, ar and esr genes, respectively. The Maximum Likelihood (ML) phylogenetic trees were generated using the IQ-Tree (http://iqtree.cibiv.univie.ac.at/ (accessed on 4 June 2021 )) with a bootstrap of 1000 replicates and were visualized by the iTOL (https://itol.embl.de/ (accessed on 6 June 2021)). The best models selected for phylogeny construction were the JTT + F + I + G4 model for cyp genes, the VT + F + G4 model for hsd genes, the JTT + G4 model for star, the JTT + F + I + G4 model for ar, and the JTT + F + I + G4 model for esr.

Molecular Evolution Analysis
Molecular evolution analysis was performed with PAML [52] in EasyCodeML V1.21 [53] to compare the evolution rates of steroidogenic genes in teleost lineages, including 10 teleost species as Japanese flounder (P. olivaceus), Chinese tongue sole (C. semilaevis), zebrafish (D. rerio), stickleback (G. aculeatus), medaka (O. latipes), fugu (T. rubripes), Asian sea bass (Lates calcarifer), platyfish (Xiphophorus maculatus), tilapia (Oreochromis niloticus) and turbot (Scophthalmus maximus). Firstly, the ratio of non-synonymous to synonymous substitutions (dN/dS, ω) and the likelihood ratio tests (LRTs) was employed to evaluate the selective pressure in each gene. Six site model (SM) tests were used to evaluate the positive selection in each codon: M0 assumes to have the same ω for all codons, M3 assumes the ω of all codons showing a simple discrete division trend; M1a assumes only conservative sites (0 < ω < 1) and neutral sites (ω = 1) for all codons, while M2a is considered to increase the existence of positive sites (ω > 1) for all codons on the basis of M1a; the ω of all codons in M7 are assumed to belong to the matrix (0, 1) with a beta distribution, while M8 adds another type of ω (ω > 1) on the basis of M7. LRTs were used to judge whether the paired models (M0 vs. M3; M1a vs. M2a; M7 vs. M8) are significantly different [52], and to estimate whether there are positively selected sites (PSSs) (M2a vs. M1a and M8 vs. M7) under the premise of a significant p-value (p < 0.05). For branch model (BM) tests, firstly, the P. olivaceus branch was labeled as the foreground branch (ω1) and the other teleosts as the background branch (ω0) to compare the evolution rates of steroidogenic genes in P. olivaceus against other teleosts; secondly, for the putative duplication originated genes, each gene cluster was labeled as the foreground branch (ω1) to compare it with the other branch (ω0) in teleost lineages.

Expression of Steroidogenic Genes through Transcriptomic Analysis
To investigate the expression profiles of steroidogenic genes, four RNA-seq datasets from P. olivaceus and C. semilaevis were obtained from NCBI. For P. olivaceus, one dataset was from the gonad tissues of 1.5-year-old wildtype female ovary (FO), wildtype male testis (MT), gynogenetic female ovary (GFO), gynogenetic neo-male testis (GNMT) (PR-JNA639001) [54], and the other dataset was from the brain and gonad tissues of 1.5-year-old gynogenetic P. olivaceus (PRJNA764760) [34]. Two transcriptomes for C. semilaevis were also employed, such as liver, brain, pituitary and gonad tissues of 2-year-old male and female individuals (PRJNA413516) [55], as well as the brain and gonad tissues of 1-year-old male, female and pseudo-male individuals (PRJNA480118) [27].
For transcriptome analysis, the fastq files were first optimized by trimmomatic [56], to remove adapters, low-quality reads, and the reads with less base numbers. Next, the read mapping and quantification were performed with the Hisat and StringTie pipeline [57], and the Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) values were extracted for each steroidogenic gene. Differentially expressed genes (DEGs) between groups were determined by edgeR [58]. Any pair-wise genes with log 2 |fold change| (log 2 FC) > 2, adjusted p-value < 0.05 and at least one gene with FPKM >1 was identified as the DEG. TBtools were utilized to draw heat maps with both log 2 (FPKM + 1) and log 2 FC values. GO and KEGG annotation of steroidogenic genes were obtained from DAVID (https://david.ncifcrf.gov/ (accessed on 8 September 2021)) for humans and zebrafish. Redundancy (semantically similar terms) and pathways involving a low number of genes (<3) were removed.

Protein-Protein Interaction and Co-Expression Analysis
Protein-protein interaction (PPI) networks of steroidogenic genes in humans and zebrafish were analyzed using the online site STRING 11.5 (https://string-db.org/ (accessed on 14 August 2021)) with the parameter, minimum required interaction score, set to 0.7. For P. olivaceus and C. semilaevis, OmicShare (https://www.omicshare.com/ (accessed on 17 August 2021)) was employed for the Pearson's correlation coefficient (PCC) analysis of FPKM values between steroidogenic genes with a PCC threshold of 0.8 from the transcriptome analysis. Cytoscape [59] was employed to illustrate the networks according to the PCC values.

Conclusions
In summary, steroidogenic genes were comprehensively characterized in P. olivaceus and C. semilaevis, with cyp19a1a/cyp19a1b, hsd17b12a/hsd17b12b, ara/arb and esr2a/esr2b originating from the TSGD, and star/star2 and cyp17a1/cyp17a2 putatively evolving independently. The expression profile and network analysis indicated the key factors functioning essentially during the steroidogenesis of P. olivaceus and C. semilaevis. Specifically, star, cyp11a/c, cyp17a1/a2, cyp21a, hsd3b1/b7, hsd11b and hsd20b were primarily expressed in the testis of P. olivaceus and C. semilaevis, whereas cyp19a1a and hsd17b genes were highly expressed in the ovaries, which constitute the coordinated regulatory network during the steroidogenic process. Interestingly, only a few steroidogenic genes were differentially expressed between male and neo-male testis, which warrants more functional verification. Investigation of other functional pathways, such as the piRNA pathway [60], may help to explain the specific difference in testis development and spermatogenesis between male and neo-male P. olivaceus and C. semilaevis in the future. All these findings will provide significant insights into the steroidogenic enzyme cascade of the teleost lineages and may contribute to reproductive manipulation in aquaculture.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ijms23105520/s1. Author Contributions: J.C. designed the study, F.Y., Y.W., W.L., Q.Z. and W.Z. analyzed the data. J.C. and F.Y. conducted the manuscript writing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The transcriptome datasets used in this study can be found in the NCBI Sequence Read Archive (SRA) BioProject as PRJNA639001 and PRJNA764760 for P. olivaceus, and PRJNA413516 and PRJNA480118 for C. semilaevis.