Two duplicated gsdf homeologs cooperatively regulate male differentiation by inhibiting cyp19a1a transcription in a hexaploid fish

Although evolutionary fates and expression patterns of duplicated genes have been extensively investigated, how duplicated genes co-regulate a biological process in polyploids remains largely unknown. Here, we identified two gsdf (gonadal somatic cell-derived factor) homeologous genes (gsdf-A and gsdf-B) in hexaploid gibel carp (Carassius gibelio), wherein each homeolog contained three highly conserved alleles. Interestingly, gsdf-A and gsdf-B transcription were mainly activated by dmrt1-A (dsx- and mab-3-related transcription factor 1) and dmrt1-B, respectively. Loss of either gsdf-A or gsdf-B alone resulted in partial male-to-female sex reversal and loss of both caused complete sex reversal, which could be rescued by a nonsteroidal aromatase inhibitor. Compensatory expression of gsdf-A and gsdf-B was observed in gsdf-B and gsdf-A mutants, respectively. Subsequently, we determined that in tissue culture cells, Gsdf-A and Gsdf-B both interacted with Ncoa5 (nuclear receptor coactivator 5) and blocked Ncoa5 interaction with Rora (retinoic acid-related orphan receptor-alpha) to repress Rora/Ncoa5-induced activation of cyp19a1a (cytochrome P450, family 19, subfamily A, polypeptide 1a). These findings illustrate that Gsdf-A and Gsdf-B can regulate male differentiation by inhibiting cyp19a1a transcription in hexaploid gibel carp and also reveal that Gsdf-A and Gsdf-B can interact with Ncoa5 to suppress cyp19a1a transcription in vitro. This study provides a typical case of cooperative mechanism of duplicated genes in polyploids and also sheds light on the conserved evolution of sex differentiation.

Introduction Polyploidy or whole-genome duplication (WGD) provides extra substrates for genomic evolution and is thus considered to be an important driving force for genetic diversity, trait innovation, and ecological adaption [1][2][3][4][5]. The majority of plants and vertebrates have evolved from polyploid ancestors [6]. Recent polyploidy is also widespread in fishes and amphibians, though it is apparently less frequent than in plants [5,7]. Polyploidy may inherit additional set/sets of chromosomes from the same species (autopolyploidy) or from interspecific hybridization (allopolyploidy) [6]. Homologous chromosomes (and the genes they carry) resulting from allopolyploidy are commonly referred to as homeologs (also homoeologs) [3]. During the initial polyploidization, the neopolyploids are thought to experience genomic chaos resulting from the emergence of duplicated genomes [8][9][10][11]. During the subsequent diploidization processes, duplicated genes generated from polyploidy will be eliminated or will evolve divergent functions, and their evolutionary fates are closely associated with the interplay of structural and functional entanglement [12,13]. Recent studies in goldfish and common carp suggest that the subgenomes resulting from the allotetraploidy have continuously rediploidized in a manner from asymmetrical evolution to diverse stabilization [14][15][16]. The homeologs from subgenomes A and B are co-expressed in most pathways, and their expression dominance shifts temporally during embryogenesis [14]. However, the molecular mechanisms underlying how the duplicated genes co-regulate a biological process in polyploids remains largely unknown.
Except for the two rounds (1R and 2R) of WGD that the common vertebrate ancestor has undergone [17], a fish-specific WGD (3R) is believed to result in the dramatic radiation of teleosts [18,19]. Intriguingly, the hexaploid gibel carp (Carassius gibelio), a cyprinid fish with a wide distribution across Eurasia [20,21], have undergone two extra rounds of polyploidy [22]. An early allotetraploidy about 10−15 Mya resulted in the formation of tetraploid C. auratus (AABB, 4n = 100) [14,23], and a late extra autotriploidy from an ancestral tetraploid approximately 0.5 Mya led to the occurrence of hexaploid gibel carp (C. gibelio) (AAABBB, 6n�150) [22,24]. The hexaploid gibel carp is actually an amphitriploid with two triploid sets of chromosomes derived from both ancestors [22,25,26]. Thus, in hexaploid gibel carp, most genes usually have a total of two homeologs, one from subgenome A and the other from subgenome B, and each homeolog commonly has three alleles, thereby providing an ideal model to investigate cooperative mechanisms of duplicated genes.
The amphitriploid gibel carp has overcome the meiotic obstacle caused by three homologous chromosomes via unisexual gynogenesis, in which the eggs are activated by the sperm of sympatric sexual species (kleptospermy) to initiate embryogenesis using only maternal genetic information [27,28]. In contrast to other unisexual vertebrates, variable male proportions ranging from 1.2% to 26.5% have been discovered in wild populations [29,30]. In our previous studies, we identified a genetic male-specific marker (MSM) in gibel carp [31] and revealed that male-specific supernumerary microchromosomes are closely associated with the occurrence of genotypic males in a dose-dependent relationship [32,33]. When a female gibel carp is mated with a male from other species, a typical gynogenesis is initiated that all the offspring have the same genetic information as the maternal individual. When a female gibel carp is mated with a genotypic amphitriploid male, a variant of gynogenesis is initiated, during which the sperm nuclei are also extruded but some supernumerary microchromosomes of sperm nuclei occasionally leak into the unreduced eggs. This variant of gynogenesis can accumulate microchromosomes, generate males, and create genetic diversity in the offspring [32][33][34][35]. However, details concerning the molecular mechanism of male determination and differentiation in this gynogenetic hexaploid fish are limited.
In sharp contrast to the remarkable diversity of sex-determining switches [36][37][38], the downstream genetic cascades of sex differentiation are relatively conserved [37,[39][40][41]. Members of the transforming growth factor-β (TGF-β) signaling pathway have been identified as being vastly involved in sex determination and differentiation in vertebrates [36,42]. Gonadal somatic cell-derived factor (gsdf), a member of the TGF-β superfamily [43], commonly acts as a male gonad factor during sex determination/differentiation in fish species [36,37,[44][45][46][47][48][49], but the molecular details underlying gsdf-mediated male determination/differentiation remain elusive. In this study, we choose the hexaploid gibel carp and the gsdf gene, as a unique system to analyze the cooperative mechanisms of duplicated genes on polyploid sex differentiation. We identified two divergent gsdf homeologous genes and revealed that the gsdf homeologs cooperatively regulated male differentiation by inhibiting cyp19a1a transcription. And interactive mechanism analyses illustrated that Gsdf-A and Gsdf-B could interact with Ncoa5 to suppress cyp19a1a transcription in vitro.

Characterization of gsdf homeologs and alleles
In the hexaploid gibel carp (C. gibelio), we identified two divergent gsdf homeologs derived from subgenome A (gsdf-A) and subgenome B (gsdf-B), which were localized to chromosomes A21 and B21, respectively ( Fig 1A). As gibel carp reproduce via unisexual gynogenesis without meiotic recombination [26,27], single nucleotide polymorphisms (SNPs) specific to each allele were stable through generations [25]. According to the SNPs, the sequenced fragments of each homeolog could be clearly divided into three types, indicating that gsdf-A and gsdf-B had three alleles each ( Fig 1B). Subsequently, we cloned the coding sequences of three gsdf-A alleles and three gsdf-B alleles, and found that the average identity between gsdf-A and gsdf-B was 85.36 ± 0.34%, while the average identities among three alleles of gsdf-A and three alleles of gsdf-B were 99.53 ± 0.27% and 99.29 ± 0.36%, respectively (S1A Fig). Interestingly, the deduced amino acid sequences of gsdf-A and gsdf-B were less conserved (average identity = 77.66% ± 0.37%) than their coding sequences (S1B Fig), as most of the differing nucleotides (68.70 ± 1.33%) caused amino acid changes (S1 Table).
Both gsdf-A and gsdf-B contain five exons and four introns, which is similar to the gsdf genes of C. auratus, Danio rerio, Ictalurus punctatus, Oryzias latipes, and Oreochromis niloticus ( Fig 1C). In addition, most genes in the neighborhood around gsdf showed conserved synteny in genomic blocks among these fish species (Fig 1D). Subsequently, we performed sequence alignment of the deduced amino acids and ascertained that the Gsdf-A and Gsdf-B of hexaploid C. gibelio had conserved TGF-β domains as in other fish species, particularly the seven or eight cysteines in this domain (S1B Fig). Phylogenetic reconstruction showed that Gsdf-A of C. gibelio and C. auratus were clustered in one clade, while Gsdf-B of C. gibelio and C. auratus were clustered into another clade (S1C Fig), and these patterns were in accordance with the common allopolyploidy origin shared by these two fish species [14,22,26].

Dynamic transcription of gsdf-A and gsdf-B are mainly activated by dmrt1-A and dmrt1-B, respectively
We first examined gsdf transcription in eight adult organs via relative real-time quantitative polymerase chain reaction (qPCR) and found out that both gsdf-A and gsdf-B mRNAs were distributed exclusively in the gonads, with much higher expression in the testis than in the ovary (Fig 2A). Then, we analyzed the dynamic expression profiles of gsdf-A and gsdf-B in male gonads during developmental stages at 17,21,30,45,60,90,120,150,210,250, 300, and 360 days post hatching (dph). In gibel carp, gonadal morphological differentiation between females and males commonly occurs around 40 dph and gonad mature at about 1 year post hatching [33,50]. The expression of gsdf-A steadily increased with a very slow growth rate before 120 dph, displayed a sharp increase from 150 to 250 dph, and then peaked at 250 dph. The expression of gsdf-B rapidly increased at 30 dph and peaked at 120 dph. After reaching the peak value, the expression levels of both gsdf-A and gsdf-B decreased in the mature testis and remained at a certain level ( Fig 2B). Subsequently, we produced a polyclonal antibody against Gsdf that could recognize both Gsdf-

PLOS GENETICS
Two gsdf homeologs regulate male differentiation by inhibiting cyp19a1a analysis was performed to assess the cellular distribution of Gsdf proteins (Gsdf-A and Gsdf-B) in male gonads at 30, 90, and 360 dph. According to the fluorescence intensity of the anti-Vasa antibody and nuclear morphology, we could easily distinguish germ cells from somatic cells [50,51]. Along with spermatogenesis, the green signal derived from Gsdf was mostly distributed in the cytoplasm of somatic cells surrounding the germ cells ( Fig 2C), which was consistent with the main distribution in Sertoli cells of other fishes [46,47,52,53]. Gsdf, as a member of the TGF-β superfamily, is usually considered to be a secreted ligand that will bind to its cell surface receptors [54,55]. However, Gsdf proteins might also accumulate in cytoplasm before being excreted.
By in silico analysis, three predicted Dmrt1-binding sites (S3 Fig) and two predicted Sf1 (also named as Nr5a1)-binding sites (S4 Fig) [56] were identified in the upstream sequences of both gsdf-A (from -2080 to +50) and gsdf-B (from -2150 to +50), which were defined as potential promoter of gsdf-A (2130 bp) and gsdf-B (2200 bp), respectively ( Fig 3A). Subsequently, we cloned these two potential promoters into a pGL3-Basic luciferase reporter vector to analyze Dmrt1's capability for activating gsdf promoters. Renilla luciferase plasmid pRL-TK was used as an internal reference. Expression plasmids of Sf1-A/Sf1-B, Dmrt1-A, and Dmrt1-B were constructed, and an empty expression plasmid was used as control. Similar to Nile tilapia [57], gibel carp gsdf transcription was also activated by Dmrt1 in a dose-dependent manner in the presence of Sf1 in Carassius auratus L. blastulae embryonic (CAB) cells ( Fig 3B). Intriguingly, for the gsdf-A potential promoter, the transcriptional regulation ability of Dmrt1-A was much stronger than that of Dmrt1-B. For the gsdf-B potential promoter, the transcriptional regulation ability of Dmrt1-B was much stronger than that of Dmrt1-A ( Fig 3B). Mutation of the Dmrt1-binding site 1 or 3 on gsdf-A potential promoter led to a decrease in both Dmrt1-Ainduced and Dmrt1-B-induced transcriptional activation of gsdf-A. Meanwhile, mutation of the Dmrt1-binding site 1 or 3 on gsdf-B potential promoter resulted in a decrease in both Dmrt1-A-induced and Dmrt1-B-induced transcriptional activation of gsdf-B. However, mutation of Dmrt1-binding site 2 on gsdf-A and gsdf-B did not significantly affect Dmrt1-A/ Dmrt1-B-induced transcriptional activation (Fig 3C and 3D). Thus, Dmrt1-binding site 1 and 3 of both gsdf-A and gsdf-B were important for gsdf activation.
Subsequently, we found out that the expression patterns of dmrt1-A and dmrt1-B in both adult organs and in gonads at different developmental stages (Fig 3E and 3F) were closely associated with the expression patterns of gsdf-A and gsdf-B, respectively (Fig 2A and 2B). Fluorescence in situ hybridization (FISH) analysis showed that dmrt1 and sf1 mRNA were coexpressed in some Gsdf positive somatic cells in testis ( Fig 3G). In addition, we performed dN/ dS analyses of dmrt1 and gsdf genes to check whether one or the other homeologs are under selection. Using D. rerio as reference, the dN/dS value of dmrt1-A was significant lower than that of dmrt1-B (χ 2 test p value: 7.93×10 −4 ) under the two-ratio model (Fig 3H and 3J). On the other side, the dN/dS analysis of gsdf homeologs met one-ratio model and display the same value (Fig 3I and 3K). These results indicate that the dmrt1 and gsdf homeologs are under asymmetric and symmetric purifying selection, respectively. Thus, the differential expression between gsdf-A and gsdf-B might be resulted from the divergent evolution of dmrt1-A and dmrt1-B.

Deficiency of gsdf-A or/and gsdf-B leads to partial/complete male-to-female sex reversal
To uncover the function of gsdf in male development, we performed loss-of-function analysis using CRISPR/Cas9 in the hexaploid C. gibelio with three alleles of gsdf-A and three alleles of   [26,27], different alleles may have different mutations in one mutant individual [25]. For instance, all the individuals in the gsdf-A mutant line (gsdf-A +1/Δ7/Δ4 + gsdf-B +/+/+ MSM+) had the same genotype where the first allele of gsdf-A had a 1-bp insertion, the second allele of gsdf-A had a 7-bp deletion, the third allele of gsdf-A had a 4-bp deletion, and three alleles of gsdf-B were all wild type genotype without mutations ( Fig 4A).

PLOS GENETICS
Two gsdf homeologs regulate male differentiation by inhibiting cyp19a1a demonstrated to play a role in female differentiation in hexaploid C. gibelio (S6 Fig). In the undifferentiated gonads of all gsdf mutants (MSM+), the ovarian aromatase gene cyp19a1a was dramatically upregulated (Fig 5C), indicating that deficiency of gsdf might abolish its repression of cyp19a1a. Thus, we used nonsteroidal aromatase inhibitor letrozole to analyze whether blockage of Cyp19a1a enzyme activity could rescue the male-to-female sex reversal caused by gsdf dysfunction. As expected, administration of letrozole from 15 to 55 dph prevented male-to-female sex reversal in 84.4% of gsdf-A/gsdf-B double mutants (MSM+) ( Fig  5F). As the period before 15 dph was key developmental stages of sex determination/differentiation, letrozole treatment after 15 dph might be the reason that 15.6% gsdf-A/gsdf-B double mutants still developed into phenotypic females. And the gene expression patterns of these letrozole treated testes were similar to that in WT males (MSM+), except for cyp19a1a and foxl2b (Fig 5G). Upregulation of cyp19a1a and its active transcriptional factor foxl2b [25,60] in the testis of gsdf-A/gsdf-B double mutants (MSM+) subjected to letrozole treatment may have been caused by inhibition of Cyp19a1a enzyme activity. In addition, cyp19a1a mRNA was coexpressed with Gsdf in some somatic cells of mature testis (S7 Fig). These results indicate that gsdf-A and gsdf-B co-inhibit cyp19a1a expression, resulting in male differentiation in WT males (MSM+).

Identification of Gsdf-Ncoa5 interaction via yeast two-hybrid assay and coimmunoprecipitation
As a member of TGF-β superfamily, Gsdf is commonly considered to be a ligand to bind to its cell surface receptors [54]. To identify potential interaction membrane proteins of secreted Gsdf, the coding sequence of Gsdf-A mature peptide was cloned into pBT3-SUC (pBT3-SUC-Gsdf-A mature peptide) as a bait to perform yeast two hybrid assay via DUAL membrane system (Dualsystems Biotech). However, screens did not yield any interactors (S8 Fig). Subsequently, the open reading frame of gsdf-A was cloned into pGBKT7 (pGBTKT7-Gsdf-A) as bait to perform yeast two hybrid assay via GAL4 system (Clontech). A total of 239 positive transformants were selected, and their plasmids were isolated for sequencing. These obtained coding sequences belonged to 34 genes according to the genome and transcription data of gibel carp [26,33]. To exclude false positive results, the full-length coding sequences of the 34 genes were cloned into vector pGADT7-AD and co-transformed with gsdf-A-pGBTKT7 separately. Finally, a total of 27 proteins were confirmed to be the potential interaction partners of Gsdf-A (Fig 6A and 6B).
We overexpressed these 27 isolated genes in order to identify whether these genes were involved in sex differentiation pathway, by analyzing transcription of sex differentiation genes such as cyp19a1a, foxl2b, dmrt1, and amh. In CAB cells, amh had no constitutive expression, so we obtained the data of the rest three sex differentiation genes (S9 Fig). Among the 27 proteins, a total of 6 proteins could activate cyp19a1a transcription significantly, however, only Ncoa5 activated cyp19a1a transcription but did not change expression levels of foxl2b and dmrt1 (S9 Fig). Besides, Ncoa5 was revealed to be involved in the regulating pathway of sex hormones [61,62] and cyp19a1a was one of the most important downstream genes of gsdf ( Fig  5), so we selected Ncoa5 for subsequent analyses. There were two homeologs of ncoa5 (ncoa5-A and ncoa5-B) in the gibel carp genome (S10A and S10C Fig) and the protein isolated from yeast two-hybrid assay was Ncoa5-B ( Fig 6B). As the protein sequences between Ncoa5-A and Ncoa5-B was highly conserved (identity = 91.53%) (S10C Fig), we only used Ncoa5-B for subsequent in vitro experiments.
Subsequently, we found that the presence of Gsdf-A/Gsdf-B could eliminate Ncoa5-induced upregulation of cyp19a1a expression (Fig 6C), but could not eliminate Foxl2b/ Sf1-induced upregulation of cyp19a1a (Fig 6D). In addition, the interactions of Gsdf-A/Gsdf-B and Ncoa5 were confirmed by co-immunoprecipitation, in which the anti-Flag Ab-immunoprecipitated protein Ncoa5 was recognized by the anti-HA Ab (Fig 6E). These findings indicate that Gsdf-A and Gsdf-B both interact with Ncoa5 to regulate cyp19a1a transcription in vitro.

Ncoa5 participates in cyp19a1a regulation via interaction with Rora
In humans, Rora is known to interact with Ncoa5 to enhance cyp19a1a transcription [61,62]. It would be interesting to know whether gibel carp Rora and Ncoa5 are involved in the expression modulation of cyp19a1a. We identified two homeologs of roraα (roraα-A and roraα-B) and two homeologs of roraβ (roraβ-A and roraβ-B) in the gibel carp genome (S10B and S10D Fig). The sequence identities between Roraα-A/Roraα-B and human RORA was much higher than those between Roraβ-A/Roraβ-B and human RORA (S10D Fig). As the protein sequences between Roraα-A and Roraα-B (identity = 91.88%) was highly conserved (S10D Fig) and Roraα-B had higher identity than Roraα-A compared with human ortholog, we used Roraα-B for subsequent in vitro experiments.
In CAB cells, co-immunoprecipitation showed that Rora could interact with Ncoa5 ( Fig  7A). In addition, Rora activated cyp19a1a transcription in a dose-dependent manner in the presence of Ncoa5 (Fig 7B). Potential Rora-binding sites were predicted in the promoter of cyp19a1a (Figs 7C and S11) [63]. Mutation of the potential Rora-binding sites resulted in a decrease in the Rora/Ncoa5-induced transcriptional activation of cyp19a1a but did not affect Foxl2b/Sf1-induced upregulation of cyp19a1a, indicating that these binding sites were specific to Rora/Ncoa5 (Fig 7D). Subsequently, the Rora-binding site 1 was further confirmed via chromatin immunoprecipitation (ChIP) in CAB cells transfected with Rora-Myc vector. The PCR band containing the binding site 1 of Rora was detected in the chromatin precipitated with the Myc antibody, while no band was observed in the negative control chromatin that was precipitated with nonspecific IgG (Fig 7E). Besides, FISH analysis showed that rora and ncoa5 mRNA were co-expressed in Cyp19a1a positive somatic cells of mature testis (Fig 7F). Thus, these findings indicate that Rora positively regulates transcription of cyp19a1a by binding to the promoter in the presence of Ncoa5.

Gsdf-A and Gsdf-B both inhibit cyp19a1a transcription via competitive interaction with Ncoa5
During sex differentiation, expression of gsdf was closely negatively associated with expression of cyp19a1a in female and male gonads (Fig 8A). The expression of ncoa5 had no significant difference between females and males in gonads at early developmental stages, while the gonadal expression of rora has no difference between females and males before 45 dph but display female-biased expression at 60 dph (Fig 8A). FISH analysis showed that rora and ncoa5 mRNA were co-expressed in some Gsdf positive somatic cells of mature testis (Fig 8B). To elucidate how gsdf-A or gsdf-B regulates cyp19a1a expression, we performed an overexpression analysis, and the overexpression of Gsdf-A or Gsdf-B significantly repressed cyp19a1a transcription in vitro accompanied by downregulation of rora (Fig 8C and 8D). However, the expression of ncoa5 was not affected by Gsdf-A or Gsdf-B overexpression in tissue culture cells (Fig 8C and 8D), consistent with the in vivo result. Meanwhile, the cyp19a1a transcription induced by Rora and Ncoa5 was repressed by Gsdf-A and Gsdf-B in a dose-dependent manner ( Fig 8E). Furthermore, siRNAs against both ncoa5-A and ncoa5-B were designed, and the repression of Gsdf-A or Gsdf-B on cyp19a1a transcription was abolished by knockdown of both ncoa5-A and ncoa5-B (Fig 8F). The competitive interaction analysis revealed that Rora recruitment by Ncoa5 decreased upon the participation of Gsdf-A or Gsdf-B but was not affected by the participation of other proteins such as female differentiation factor Foxl2b and immune-related factor Sting ( Fig 8G). Thus, these results indicate that the upregulation of Gsdf-A and Gsdf-B increases the interaction of Ncoa5 with both homeologs and reduces the interaction of Ncoa5 with Rora, leading to the downregulation of cyp19a1a and subsequent male differentiation.

Discussion
Polyploidy provides a source of new genes and these duplicated genes will be eliminated/pseudogenized or evolve a sub/neo function during the evolutionary trajectory [2,12,25]. The

PLOS GENETICS
Two gsdf homeologs regulate male differentiation by inhibiting cyp19a1a hexaploid gibel carp (AAABBB) with extra two rounds of polyploidy origins has retained most of the duplicates, where most genes usually have two homeologs, and each homeolog commonly has three alleles [22,25,26,64]. Recently, asymmetrical evolution, homoeologous exchanges, and expression divergence of subgenomes A and B have been observed in allotetroploid goldfish, common carp, and hexaploid gibel carp [14,65]. And we also have presented functional divergence of foxl2 and viperin homeologs in hexaploid gibel carp [25,64]. Here, we found expression divergence between gsdf-A and gsdf-B (Fig 2A and 2B), but knockout gsdf-A or gsdf-B displayed similar sex-reversal rates (S5 Fig), indicating the contribution of each homeolog is similar in male differentiation. In addition, disruption of gsdf-A or gsdf-B triggers highly compensatory expression of gsdf-B or gsdf-A during the critical period of sex differentiation (Fig 5E), and missing either gsdf-A or gsdf-B does not give complete sex reversal but missing both does, suggesting that gsdf-A and gsdf-B cooperatively regulate male differentiation in gibel carp.
An intriguing finding of this study is the revelation of potential molecular rationales underlying male differentiation mediated by gsdf homeologs. In male gibel carp with MSM (MSM+), high expression of gsdf-A and gsdf-B in somatic cells suppresses cyp19a1a to induce Sertoli cell development and male differentiation. In female individuals without MSM (MSM−), the low levels of Gsdf-A and Gsdf-B cannot inhibit cyp19a1a transcription, leading to estrogen production, granulosa cell development, and female differentiation (Fig 9). Besides, in vitro analyses revealed that Gsdf-A and Gsdf-B can interact with Ncoa5 to block Ncoa5 interaction with Rora, inhibiting Rora/Ncoa5-induced activation of cyp19a1a (Fig 9). Commonly, as a member of TGF-β superfamily, Gsdf is considered to be a secreted ligand to bind to its cell surface receptors [54]. Here we demonstrated that Gsdf might also have functions in cells. However, we still do not known whether other factors are involved in the interaction between Gsdf and Ncoa5. As we known, Foxl2/Sf1-cyp19a1a pathway is an important pathway of female sex differentiation and knockout of foxl2 results in female to male sex reversal in fish [25,66]. We suppose that the pathway of Rora/ Ncoa5 induced cyp19a1a activation is not independent of Foxl2/Sf1-cyp19a1a pathway. For instance, the cyp19a1a expression change induced by Foxl2/Sf1 can affect estradiol synthesis, which may lead to expression changes of er (estrogen receptor) and rora in the presence of Ncoa5 and then affect cyp19a1a transcription via Rora/Ncoa5-cyp19a1a pathway [61,62]. Thus, these two pathways of cyp19a1a regulation may interact with each other.
The sexual phenotype is the result of antagonism between the female and male pathways, with multiple feedback loops that are influenced by genotypic and/or environmental factors [37]. In hexaploid gibel carp, we have found that gsdf and cyp19a1a play antagonistic roles in sex differentiation. Gsdf represses cyp19a1a by blocking Ncoa5's availability for activation of cyp19a1a transcription. Meanwhile, cyp19a1a also has the ability to downregulate gsdf by suppressing gsdf's transcription factor dmrt1 [59]. Antagonistic actions of dmrt1 and foxl2 have been found in many other vertebrates [67][68][69][70] and these two genes also display conserved expression patterns during sex differentiation in gibel carp (Fig 5) [25,71]. In addition, Foxl2 also has been demonstrated to be a positive transcriptional factor of cyp19a1a as previously reported [60], indicating that Dmrt1 and Foxl2 may also play conserved antagonistic roles in gibel carp. As shown in previous studies, on one hand upregulation of cyp19a1a would lead to dmrt1 inhibition while on the other it would upregulate foxl2.
In this study, we have identified two duplicated gsdf homeologous genes, gsdf-A and gsdf-B, and each homeolog has three alleles in the gynogenetic hexaploid gibel carp. The transcription of gsdf-A and gsdf-B is mainly activated by dmrt1-A and dmrt1-B, respectively. Moreover, lossof-function experiments reveal the cooperative ability of two gsdf homeologs to regulate male differentiation by interacting cyp19a1a transcription. And the interactive mechanism analyses demonstrate that Gsdf interacts with Ncoa5 to suppress cyp19a1a transcription in vitro. This study provides a typical case of cooperative mechanism of duplicated genes in polyploids and also sheds light on the conserved evolution of sex differentiation.

Ethics statement
Animal experiments and treatments were performed according to the Guidelines for Animal Care and Use Committee of Institute of Hydrobiology, Chinese Academy of Sciences (IHB, CAS, Protocol No. 2016-018).

Fishes and cells
Experimental fish species including hexaploid gibel carp (C. gibelio) and red common carp (C. carpio) were provided and raised by the National Aquatic Biological Resource Center (NABRC), Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China. Fish cell line Carassius auratus L. blastulae embryonic (CAB) cells and epithelioma papulosum cyprini (EPC) cells were maintained at 28˚C in 5% CO2 in medium 199 (Invitrogen) supplemented with 10% fetal bovine serum (FBS) (Invitrogen).

Cloning and sequence analysis
The divergent gsdf homeologs including gsdf-A and gsdf-B were identified according to the assembled genome of hexaploid gibel carp (C. gibelio) (Genbank accession numbers: PRJNA546443). Full-length cDNAs of gsdf-A and gsdf-B were obtained by 5' and 3' rapid amplification of cDNA ends (RACE) (SMARTer RACE 5'/3' Kit, Clontech) using testicular cDNA library. Specific primers for RACE amplification (S2 Table) were designed according to the genome sequences of gsdf-A and gsdf-B. Multiple alignments of gsdf genomic and cDNA sequences was performed by DNAMAN 8.0 software.
The deduced amino acid sequences were predicted by DNAMAN 8.0 software. Multiple alignments of deduced amino acid sequences was performed by ClustalX program and exhibited by Bioedit program. Phylogenetic construction was adjusted by bootstrap analysis (1000 replicates) using the neighbor-joining method (NJ) in MEGA version 7.0 [72].

RNA extraction and qPCR
Adult organs, including heart, liver, hypothalamus, pituitary, kidney, spleen, ovary and testis, were isolated from three mature WT females (MSM-) and three mature WT males (MSM+), respectively. Male individuals were obtained from the offspring of a WT female mating with a WT male and all the female offspring were excluded by PCR detection using the male-specific marker. WT male gonads at different developmental stages were carefully dissected, and a total of 25, 25, 20, 10, 5, 5, 5, 5, 5, 5, 5, and 5 gonads were pooled for RNA extraction at the developmental stages of 17, 21, 30, 45, 60, 90, 120, 150, 210, 250, 300, and 360 dph, respectively. Besides, a total of 25, 5, 5 gonads from gsdf mutants and corresponding WT individuals were pooled for RNA extraction at the stages of 25, 55, 100 dph, respectively.
Total RNA isolation was performed using SV Total RNA isolation System (Promega), and the isolated RNAs were reverse-transcribed by the PrimeScript RT Reagent Kit (Takara). qPCR was performed on S1000 Thermal Cycler (BioRad), using iQSYBR Green Supermix (BioRad) as described previously [73]. β-actin was used as internal reference. All samples were analyzed in triplicates, and relative expression level of target gene was calculated with 2 -ΔΔCT method. The highest expression level in each qPCR analysis was used as control and defined as 1 separately. Data were displayed as mean ± standard deviation. Significant differences were calculated by one-way ANOVA followed by Tukey test.

Polyclonal anti-Gsdf antibody preparation and western blot
The gsdf-A cDNA sequence coding for 170 amino acids (S1 Fig) was cloned into Frd-GST vector (Friendbio Science and Technology) and the prokaryotic fusion protein was used as antigen to immunize a rabbit. Anti-Gsdf polyclonal antibody was produced by Friendbio Science and Technology Company Limited (Wuhan). Sample protein was extracted from cells using RIPA Lysis Buffer (Beyotime). Western blot detection was performed according to the previous reports and β-actin was used as internal control [74]. The images were obtained by Image-Quant LAS 4000mini (GE).

Histological analysis and immunofluorescence
The gonads of gibel carp were fixed with 4% paraformaldehyde in PBS at 4˚C over night. After washing with PBS, the samples were immersed in 30% saccharose-PBS buffer for 5 h at 4˚C, embedded in paraffin, and then were cut into 4μm sections. Hematoxylin-eosin staining was performed as described previously [50]. Immunofluorescence co-localization of Gsdf and Vasa was performed as described previously [71]. The images were obtained by upright fluorescence microscope Axio Imager M2 (Carl Zeiss).

Fluorescence in situ hybridization (FISH)
Probes for ncoa5 and dmrt1 antisense/sense digoxigenin-labeled RNA strands were transcribed in vitro using the DIG RNA labeling kit (Roche). Probes for roar, cyp19a1a, and sf-1 antisense/sense fluorescein-labeled RNA strands were transcribed in vitro using the Fluorescein RNA labeling kit (Roche). Specific primers with a T7 RNA polymerase promoter were designed to amplify complementary DNA (cDNA) fragment of each gene (S2 Table). Each probe was used at a final concentration of 0.5ng/μL. For more sensitive fluorescence in situ hybridization detection, the tyramide signal amplification TSA Plus Cyanine 3/Cyanine 5 System (PerkinElmer Life Science) was used according to the manufacturer's instructions. Digoxigenin-labeled RNA was stained with cy5, fluorescein -labeled RNA was stained with cy3. FISH analyses using sense RNA strands were shown in S12 Fig.

Evolutionary analysis of gsdf or dmrt1 homeologs
To investigate the potential role of selection on the evolution of gsdf or dmrt1 homeologs, the gsdf or dmrt1 gene dataset was assessed by branch model tests [75]. Alternative branch models, which allow foreground and background lineages to evolve differently (with different dN/dS), were compared to null models that assume the same ratio for all branches. Tree with divergence time was taken from MEGA version 7.0 [76]. The alternative models were evaluated for statistically significance (P < 0.05) by likelihood ratio tests (LRTs), with the null model using a χ2 distribution [77].

Generation of gsdf mutants by CRISPR/Cas9
Mutant line of gsdf-A and gsdf-B were generated by CRISPR/Cas9 as described previously [25] and the process was shown in S5 Fig. The sgRNA target sites of gsdf-A and gsdf-B were designed on the first exon and the second exon, respectively (Fig 1A and 1B). gRNAs were transcribed with the TranscriptAid T7 High-Yield Transcription Kit (Thermo Fisher Scientific). The gRNA and Cas9 protein (Invitrogen) were co-injected into one-cell-stage embryos at a concentration of 200 ng/μL and 100 ng/μL, respectively.

DNA extraction and PCR detection of MSM
A small piece of fin was used to extract genomic DNA for each sampled fish, using DNA extraction kit (Promega) according to the manufacturer's instructions. The MSM was detected by PCR using the primer pair Cg-MSM-F and Cg-MSM-R [32]. PCR analysis was performed as previously described [35].

Aromatase inhibitor treatment
Individuals from a gynogenetic family (WT, MSM−) and individuals from a gsdf-A/gsdf-B double mutant family (gsdf-A Δ4,+9/Δ4/Δ7 + gsdf-B chimeric mutations , MSM+) were divided into two groups respectively, including a control group and a treatment group. The treated fish fry were fed with fairy shrimp that had been placed in 95% ethanol containing Letrozole (MCE) at a final concentration of 150 μg/L for 0.5 hr, whereas the control fish were fed with fairy shrimp that had been socked with 95% ethanol only. The treatment lasted for 40 days from 15 to 55 dph, and then all groups were fed with normal diet and maintained in outdoor tanks as described previously [32,78].

Transient transfection
CAB cells were cultured in 6-well plates of phenol red-free M199 media (Gibco) supplemented with 10% charcoal dextran-treated serum (BI) until the cultures became approximately 75% confluent. Confluent cells were transfected using FuGENE HD Transfection Reagent (Promega) with 2 μg expression vectors. To analyze how Gsdf inhibits cyp19a1a transcription, we added 17β-estradiol (E2; Sigma-Aldrich) to the cells at a final concentration of 10 nM at 4 h post-transfection to elevate cyp19a1a transcription. The cells were harvested for RNA extraction and subsequent qPCR analysis. All the experiments were performed in triplicates.

Luciferase activity assays
CAB cells (or EPC cells) were seeded in 24-well plates of phenol red M199 medium and cotransfected with various plasmid constructs at a ratio of 10:10:1 (250 ng luciferase reporter gene plasmid: 250 ng expression plasmid: 25 ng Renilla luciferase plasmid pRL-TK) using FuGENE HD Transfection Reagent (Promega). Then, transfected cells were harvested at 24 h post-transfection and measured by the Dual-Luciferase Reporter Assay System (Promega). Luciferase activities were measured by a Junior LB9509 luminometer (Berthold, Pforzheim, Germany) and normalized to the amounts of Renilla luciferase activities. All experiments were performed at least 3 times and the significant differences were calculated by SPSS soft-ware (SPSS Inc.).
Then, these cDNAs were transferred from plasmid pDONR222 to pPR3-N-DEST (Dualsystems Biotech) by BP Clonase II enzyme mix (Invitrogen) and the library titre was 3.40 × 10 7 cfu (colony-forming units) (S13A, S13C and S13F Fig). The mature peptide sequence of gsdf-A (from 95 to 186 amino acids) was cloned into pBT3-SUC vector (Dualsystems Biotech) as bait (pBT3-SUC -Gsdf-A) and yeast library screening was performed according to protocol of the DUAL membrane starter kits User Manual (Dualsystems Biotech).
Yeast two-hybrid assay using GAL4 system cDNAs were transferred from plasmid pDONR222 to pGADT7-DEST via homologous recombination by LR Clonase II Mix (Invitrogen) (S13A, S13D and S13G Fig). The extracted pGADT7-DEST-cDNA plasmids were transfected into yeast competent cell Y187, the library titre was 1.60 × 10 7 cfu. The library titer was calculated as described previously. Insert size was identify by PCR using primer pair pGADT7-F (T7) and pGADT7-R (ADR) (S2 Table). Subsequently, the well-distributed yeast cells were cultured on 100 plates (150 mm) with dropout medium (SD/-Leu) at 30˚C for 5 days. All the appeared colonies were collected into freezing medium (YPDA medium with 25% glycerol) and stored at −80˚C for yeast two-hybrid screening. The above construction of cDNA libraries were performed by OE BioTech (Shanghai, China).
The coding sequence of gsdf-A was cloned into pGBKT7 vector (Clontech) as bait (pGBTKT7-Gsdf-A), and yeast library screening was performed according to protocol of the Matchmaker Gold Yeast Two-Hybrid System (Clontech). All positive colonies were collected from quadruple dropout medium (QDO supplemented with X-alpha-Gal and Aureobasidin A, lack of Ade, His, Leu, and Trp) separately and used for plasmid extraction. Each plasmid was transformed into TOP10 chemically competent cell for subsequent sequencing. To exclude false positive results, the full-length coding sequences of all candidate genes were constructed into pGADT7 (Clontech) and then co-transformed to Y2H competent cell with pGBKT7-Gsdf-A on QDO/X-alpha-Gal/AbA plates and DDO (lack of Leu and Trp) plates, respectively.

RNA interference
EPC cells were cultured in 12-well plates overnight, and then transfected with 50 nM small interfering RNAs (siRNA) of ncoa5 (ncoa5-A and ncoa5-B) and the negative control (si-Nc) by using FishTrans (MeiSenTe Biotechnology). 17β-estradiol was added to the cells at a final concentration of 10 nM at 4 h post-transfection to elevate cyp19a1a transcription. siRNA of ncoa5 (ncoa5-A and ncoa5-B) and si-Nc were synthesized by Sangon Biotech (Shanghai). The following sequences were targeted for ncoa5 (ncoa5-A and ncoa5-B): si-ncoa5: CCGUCAUAGUC GUCAACAATT.

Co-immunoprecipitation assay
CAB cells (or EPC cells) were seeded in 10 cm 2 dishes overnight and then transfected with a total of 10 μg of various plasmid combinations. The transfected cells were washed twice with 10 mL ice-cold PBS and then lysed by radioimmunoprecipitation (RIPA) lysis buffer with protease inhibitor cocktail (Sigma-Aldrich). After removing cellular debris, the supernatant was transferred to a 1.5 mL clean tube and incubated with 25 μL anti-Flag Affinity gel (Sigma-Aldrich) or anti-HA Magnetic Beads (Thermo Fisher) overnight at 4˚C with constant rotating incubation. Immunoprecipitated proteins were collected by Magnetic Stand (Promega), washed five times with lysis buffer, and resuspended in 100 μL SDS-PAGE protein loading buffer (Beyotime). The immunoprecipitates and whole cell lysates (WCLs) were separated by 10-12% SDS-PAGE and then transferred to polyvinylidene fluoride membranes (Millipore) for subsequent western blot analysis. Antibodies were diluted as follows: anti-β-actin (Cell Signaling Technology) at 1:3,000, anti-Flag/HA antibody (Cell Signaling Technology) at 1:3,000, anti-Myc antibody (Abcam) at 1:2,000, and HRP-conjugated anti-rabbit IgG (Thermo Scientific) at 1:5,000. Images were captured by ImageQuant LAS 4000mini (GE). Results were representative of three independent experiments.

Chromatin immunoprecipitation (ChIP)
CAB cells were cultured in 15 cm 2 plate overnight and then transfected with 20 μg expression vector pCS2+-Rora-Myc. The transfected cells were used for ChIP analysis by ChIP-IT Express Chromatin Immunoprecipitation Kits (Active Motif). The protein and chromatin complexes were immunoprecipitated by 3 μg Anti-Myc antibody (Abcam) or anti-IgG antibody (Dia-An Biotec), respectively. The purified DNA was used for PCR analysis after Immunoprecipitation purification. The PCR products were electrophoresed in 2% agarose gels. The primer pair Chip-rora-F and Chip-rora-R (S2 Table) was used to amplify specific region spanning the potential binding site for Rora.