Emilie Dujardin
email: emilie.dujardin@inrae.fr

Marjolaine André

Aurélie Dewaele

Béatrice Mandon-Pépin

Francis Poulat

Anne Aubert-Frambourg

Dominique Thepot Thépot

Investigation, Meth- Luc Jouneau

Geneviève Jolivet

Eric Pailhoux
email: eric.pailhoux@inrae.fr

Anne Frambourg

Maëlle Pannetier

DMRT1 is a testis-determining gene in rabbits and is also essential for female fertility

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

DMRT1 (Doublesex and Mab-3 Related Transcription factor 1) belongs to the highly conserved family of DM domain proteins, which exhibits a zinc finger DNA-binding motif that was initially identified in Drosophilia and Caenorhabditis elegans [START_REF] Erdman | The Drosophila doublesex proteins share a novel zinc finger related DNA binding domain[END_REF][START_REF] Raymond | Evidence for evolutionary conservation of sex-determining genes[END_REF]. Some of its orthologs have been described as Testis-Determining Factor (TDF) in vertebrate species such as medaka (Oryzias latipes) [START_REF] Matsuda | DMY is a Y-specific DM-domain gene required for male development in the medaka fish[END_REF], xenope (Xenopus laevis) [START_REF] Yoshimoto | Opposite roles of DMRT1 and its W-linked paralogue, DM-W, in sexual dimorphism of Xenopus laevis: implications of a ZZ/ZW-type sex-determining system[END_REF], or chicken [START_REF] Smith | The avian Z-linked gene DMRT1 is required for male sex determination in the chicken[END_REF]. In the last, the Z chromosome carries the DMRT1 gene. In ZZ males, two copies of the DMRT1 gene are required to induce testis determination. In ZW females and ZZ chickens harboring a non-functional copy, gonads differentiate as ovaries showing that sex determination is based on DMRT1 dosage [START_REF] Ioannidis | Primary sex determination in birds depends on DMRT1 dosage, but gonadal sex does not determine adult secondary sex characteristics[END_REF].

In mammals, where the sex-determining system is XX/XY, the TDF is the SRY gene (Sexdetermining Region of the Y chromosome) carried by the Y chromosome. Based on the mouse species, DMRT1 does not appear to have retained a crucial function in testis determination since targeted deletion of Dmrt1 only affects post-natal testis function. In fact, DMRT1 has roles in both germ cells and supporting cells in the testis, and Dmrt1 -/-males showed spermatogenesis failure with spermatogonia that did not undergo meiosis [START_REF] Matson | The mammalian doublesex homolog DMRT1 is a transcriptional gatekeeper that controls the mitosis versus meiosis decision in male germ cells[END_REF]. However, specific knock-out of Dmrt1 in adult Sertoli cells led to their transdifferentiation into granulosa cells [START_REF] Matson | DMRT1 prevents female reprogramming in the postnatal mammalian testis[END_REF]. Although DMRT1 is not required for testis determination in mice, it retained part of its function in adulthood when it is necessary to maintain Sertoli cell identity. In ovarian differentiation, FOXL2 (Forkhead family box L2) showed a similar function discrepancy between mice and goats as DMRT1 in the testis pathway. In the mouse, Foxl2 is expressed in female-supporting cells early in development but does not appear necessary for fetal ovary differentiation [START_REF] Uda | Foxl2 disruption causes mouse ovarian failure by pervasive blockage of follicle development[END_REF]. On the contrary, it is required in adult granulosa cells to maintain female-supporting cell identity [START_REF] Ottolenghi | Foxl2 is required for commitment to ovary differentiation[END_REF][START_REF] Uhlenhaut | Somatic sex reprogramming of adult ovaries to testes by FOXL2 ablation[END_REF]. In other mammalian species, such as goats, FOXL2 was shown to be crucial for ovarian determination. Indeed, naturally observed in the PIS (Polled Intersex Syndrome) mutation [START_REF] Pailhoux | A 11.7-kb deletion triggers intersexuality and polledness in goats[END_REF] or experimentally induced by genome editing in goats [START_REF] Boulanger | FOXL2 Is a Female Sex-Determining Gene in the Goat[END_REF], FOXL2 loss-of-function led to female-to-male sex reversal with the early development of XX testes. Following FOXL2 absence of expression in the XX mutant gonads (XX PIS -/-or XX FOXL2 -/-), DMRT1 was up-regulated within days before increased SOX9 expression, which then directs the differentiation of Sertoli cells and the formation of testicular cords [START_REF] Elzaiat | High-throughput sequencing analyses of XX genital ridges lacking FOXL2 reveal DMRT1 up-regulation before SOX9 expression during the sex-reversal process in goats[END_REF]. These observations in the goat suggested that DMRT1 could retain function in SOX9 activation and, thus, in testis determination in several mammals. In humans, a few eLife digest Animals that reproduce sexually have organs called gonads, the ovaries and testes, which produce eggs and sperm. These organs, which are different in males and females, originate from the same cells during the development of the embryo. As a general rule, the chromosomal sex of an embryo, which gets determined at fertilization, leads to the activation and repression of specific genes. This in turn, controls whether the cells that will form the gonads will differentiate to develop testes or ovaries.

Disruption of the key genes involved in the differentiation of the gonads can lead to fertility problems, and in some cases, it can cause the gonads to develop in the 'opposite' direction, resulting in a sex reversal. Identifying these genes is therefore essential to know how to maintain or restore fertility.

DMRT1 is a gene that drives the differentiation of gonadal cells into the testicular pathway in several species of animals with backbones, including species of fish, frogs and birds. However, its role in mammals -where testis differentiation is driven by a different gene called SRY -is not well understood. Indeed, when DMRT1 is disrupted in male humans it leads to disorders of sex development, while disrupting this gene in male mice causes infertility. To obtain more information about the roles of DMRT1 in mammalian species, Dujardin et al. disrupted the gene in a third species of mammal: the rabbit.

Dujardin et al. observed that chromosomally-male rabbits lacking DMRT1 developed ovaries instead of testes, showing that in rabbits, both SRY and DMRT1 are both required to produce testes. Additionally, this effect is similar to what is seen in humans, suggesting that rabbits may be a better model for human gonadal differentiation than mice are. Additionally, Dujardin et al. were also able to show that in female rabbits, lack of DMRT1 led to infertility, an effect that had not been previously described in other species.

The results of Dujardin et al. may lead to better models for gonadal development in humans, involving DMRT1 in the differentiation of testes. Interestingly, they also suggest the possibility that mutations in this gene may be responsible for some cases of infertility in women. Overall, these findings indicate that DMRT1 is a key fertility gene. [START_REF] Chauhan | Novel Heterozygous Genetic Variants in Patients with 46,XY Gonadal Dysgenesis[END_REF][START_REF] Ledig | Partial deletion of DMRT1 causes 46,XY ovotesticular disorder of sexual development[END_REF][START_REF] Mello | Novel DMRT1 3'UTR+11insT mutation associated to XY partial gonadal dysgenesis[END_REF]. In particular, a heterozygous de novo point mutation in the DMRT1 gene has been identified in a 46, XY individual with complete gonadal dysgenesis [START_REF] Murphy | An ancient protein-DNA interaction underlying metazoan sex determination[END_REF], suggesting that DMRT1 and SRY may be involved in testicular determination.

To clarify DMRT1 functions in non-rodent mammals, we have chosen the rabbit model, where we generated a DMRT1 mutant line thanks to the CRISPR/Cas9 technology. Firstly, we characterized the DMRT1 expression in control gonads, showing that both XY and XX fetal gonads were expressing DMRT1 before their sexual differentiation. In XY fetuses, DMRT1 and SRY presented partially overlapping territory, and somatic cells expressing both of them harbored SOX9 expression and differentiated into Sertoli cells. Secondly, thanks to our CRISPR/Cas9 genetically modified rabbit model, we demonstrated that DMRT1 was required for testis differentiation since XY DMRT1 -/-rabbits showed early male-to-female sex reversal with differentiating ovaries and complete female genitalia. However, germ cells failed to undergo meiosis, and follicles did not form in XY and XX DMRT1 -/-mutant ovaries, leading to female infertility. Finally, we demonstrated that DMRT1 was a testis-determining factor in mammals and that it was also required for female fertility.

Results

DMRT1 is expressed in genital crests of both sexes and just after SRY in XY developing testes DMRT1 expression pattern has already been reported by molecular analysis in the rabbit species from 14 days post-coïtum (dpc) to adulthood [START_REF] Daniel-Carlier | Gonad differentiation in the rabbit: evidence of species-specific features[END_REF]. We aimed to investigate further the location of the DMRT1 expression during gonadal development, firstly at earlier stages of genital crest formation (12-13 dpc; Figure 1A) using in situ hybridization (ISH). SRY expression was already detected at 12 dpc and, as expected, was found only in the XY genital ridges, where it was restricted to the medullary part of the gonad (Figure 1B). In contrast, DMRT1 was faintly expressed in the gonads of both sexes, in a few cells of the medulla under the coelomic epithelium (Figure 1B). At 12 dpc, only very few germ cells, expressing POU5F1, have completed their migration into the genital ridges (Figure 1B). Twenty-four hours later, at 13 dpc, the genital ridges had tripled in size in both sexes, and the territory of SRY expression increased within the XY developing testes (Figure 1C). The number of somatic cells expressing DMRT1 was also strongly increased in both sexes, with few of them located in the coelomic epithelium (Figure 1C). In addition, more POU5F1-expressing germ cells were detected (5-12 per section instead of 1 or 2 at 12 dpc) (Figure 1B,C).

SOX9 is detected in XY medullar cells co-expressing SRY and DMRT1

At 14 dpc, the SOX9 protein was immunodetected in a few cells located in the medullary part of the XY gonad (Figure 2). Numerous somatic cells of this region also expressed SRY and DMRT1 (Figure 2), and a few co-expressed SOX9 and DMRT1 simultaneously (Figure 2-figure supplement 1). In contrast, coelomic epithelial cells only expressed DMRT1 (Figure 2).

At 15 dpc, Sertoli cells that co-express SRY, DMRT1, and SOX9 began to be organized into embryonic cords (Figure 2). At this stage, coelomic epithelial cells expressed DMRT1 but were negative for SRY and SOX9. Furthermore, we observed an islet of cells expressing SRY and DMRT1 located in the mesonephros below the boundary with the gonad (Figure 2, dotted line). These cells expressed PAX8 (Figure 2-figure supplement 2) and could correspond to the recently described supporting-like cell population contributing to the rete testis in mice [START_REF] Mayère | Origin, specification and differentiation of a rare supporting-like lineage in the developing mouse gonad[END_REF]. As in mice, these cells will express SOX9 at the latter stages (a few of them are already SOX9 positive at 15 dpc), but unlike mice, they express SRY.

From 16 to 18 dpc, the development of the testicular cords proceeded. At these two stages (16 and 18 dpc), SRY, DMRT1, and SOX9 were expressed only in the Sertoli cells, where SRY expression began to decrease from 18 dpc (Figure 2). No more DMRT1 expression could be seen in the coelomic epithelial cells, but the tunica albuginea begins to form (Figure 2), and consequently, the coelomic epithelium will become the surface epithelium.

Persistent expression of DMRT1 in XX gonadal somatic cells until ovigerous nest formation

As described above, DMRT1 expression started at 12 dpc in the gonadal somatic compartment of both sexes (Figure 1B). In the female gonads, DMRT1 remained expressed in all somatic cells, including those of the coelomic epithelium, until 16 dpc (Figure 3A). Interestingly, as in XY gonads, we observed PAX8-positive cells in XX gonads at 15 dpc (Figure 2-figure supplement 1). These cells could contribute to the formation of the rete ovarii as in mice [START_REF] Mayère | Origin, specification and differentiation of a rare supporting-like lineage in the developing mouse gonad[END_REF]. At 18 dpc, DMRT1 expression decreased but persisted in some cells located in the coelomic epithelium and just below it, where ovigerous nest formation occurred. Interestingly, the female DMRT1-antagonist gene FOXL2 began to be expressed between 16 and 18 dpc when DMRT1 expression decreased (Figure 3B). Thereafter, at 20 dpc, DMRT1 expression was limited in some somatic cells enclosed in nascent ovigerous nests where some germinal cells also began to be positive for DMRT1 (Figure 3C and Figure 3-figure supplement 1). At this stage, DMRT1-positive territory seems to overlap that of RSPO1 but not that of FOXL2 located in the loose conjunctive tissue around the ovigerous nests (Figure 3C).

The testicular formation is impaired in DMRT1 knock-out XY rabbits

To determine the role of DMRT1 in the rabbit species used as a non-rodent mammalian model, we engineered a DMRT1 knock-out line using the CRISPR/Cas9 system with two RNA guides located in exon 3. The mutation carried by this line is a 47-bp duplication in sense, leading to a frameshift of the open reading frame and a premature stop codon (Figure 4-figure supplement 1A). This mutation does not affect DMRT1 transcription but induces a total absence of protein as shown in post-natal gonads by western blot (Figure 4-figure supplement 1B, C). Thanks to this line, we first analyzed gonadal formation at 20 dpc, when the testis and ovary were distinguishable in control animals. Indeed, at this stage, testes appeared with well-formed seminiferous cords, and ovigerous nest formation was clearly in progress in the ovaries (Figure 4A). At 20 dpc, XY DMRT1 -/-gonads failed to engage testicular differentiation and appeared quite like control ovaries, but ovarian differentiation did not appear to be affected by the loss of DMRT1 (Figure 4A). To better characterize the DMRT1 -/-gonads in XY and XX fetuses, we established the gonadal transcriptome by RNAsequencing. Heatmap representation of the 3640 differentially expressed genes in at least one of the four genotypes (adjusted p-value <0.05 and |log2FC| > 1; Supplementary file 1) was clustered into eight groups (#1 to #8, Figure 4B and Supplementary file 2). Clusters #1 and #7 contained 1331 and 315 genes, respectively, which were preferentially expressed in XY control testes. Expression of these genes was decreased in XY DMRT1 -/-gonads, harboring levels close to that of the female's ovaries (XX control or DMRT1 -/-). On the other hand, clusters #2, #3, and #5 (537, 582, and 464 genes, respectively) were composed of genes preferentially expressed in XX control ovaries, and their expression was increased in XY DMRT1 -/-gonads. Deep-sequencing transcriptomics confirmed the ovarian fate of XY DMRT1 -/-gonads. The heatmap in Figure 4C also illustrates the expression for selecting some of the main genes involved in sex determination (Figure 4C). Expression levels and patterns of the principal actors of gonadal differentiation were confirmed by quantitative RT-PCR (RT-qPCR), and the location of positive cells was achieved by immunohistochemistry. As expected, SOX9, AMH, and DHH expression levels were decreased in XY DMRT1 -/-gonads, remaining like those detected in control or DMRT1 -/-XX ovaries, while SRY expression was enhanced in XY DMRT1 -/-gonads (Figure 5A). Interestingly, we noticed a slight increase of SOX9-positive cells in XY DMRT1 -/-gonads compared to XX control or mutant ovaries (Figure 5C). In contrast, FOXL2 and CYP19A1 expression were increased in XY DMRT1 -/-gonads to similar levels to those detected in control or mutant ovaries (Figure 5B). By immunohistochemistry, we detected cells expressing FOXL2 in XY DMRT1 -/-gonads (Figure 5C). Moreover, RSPO1 expression was increased in XY DMRT1 -/- gonads, but it remained lower than in control ovaries or in XX DMRT1 -/-gonads. In the latter, the RSPO1 expression was also lower than in control ovaries, suggesting a regulatory link between DMRT1 and RSPO1 in the female pathway (Figure 5B).

Germ cells failed to engage meiosis in DMRT1 mutant gonads

After the sex determination process and the first stages of gonad formation, DMRT1 -/-gonads engage a female fate and differentiate as ovaries, whatever their sex-chromosome constitution, XX or XY. Whereas the DMRT1 expression began at 18 dpc in the XY germinal lineage of control gonads and 20 dpc in XX (Figure 3-figure supplement 1), its expression was abolished in both somatic and germ cells in DMRT1 -/-mutant gonads (Figure 5C). Although XX or XY DMRT1 -/-gonads continue to develop as ovaries, most germ cells did not engage in the meiotic process. Indeed, in control ovaries at 3 days post-partum (dpp), most germ cells were in the zygotene stage, showing nuclei with highly condensed chromatin [START_REF] Daniel-Carlier | Gonad differentiation in the rabbit: evidence of species-specific features[END_REF]; Figure 6) and were positives for Ki67, showing their exit from the G0 phase of the cell cycle (Figure 6-figure supplement 1). In contrast, in DMRT1 -/-gonads, few germ cells in the preleptotene stage were observed (Figure 6), and the majority did not express Ki67 but continued to express the pluripotency marker POU5F1 (Figure 6figure supplement 1). Subsequently, the rupture of ovarian nests and the follicle formation did not occur in DMRT1 -/-gonads. At 18 dpp, folliculogenesis had already started in control ovaries, where the first primordial follicles were visible in the deepest cortical part close to the medulla (Figure 6). In contrast, DMRT1 -/-gonads seemed to be blocked at a pre-meiotic stage, and folliculogenesis failed to occur (Figure 6). In adults, DMRT1 -/-gonads were reduced in size (Figure 6-figure supplement 2), no germ cells were detected, and some somatic cells evolved toward luteinized cells (Figure 6). Consequently, both XY and XX females were completely infertile in adulthood.

Discussion

Our study gave new insights into the conservation of the sex-determination genetic cascade across evolution. Although the signal controlling this process could take different forms in metazoans, several downstream transcription factors involved in gonadal differentiation have been conserved throughout evolution. For instance, SOX9, well known in vertebrates as being essential for Sertoli cell differentiation [START_REF] Chaboissier | Functional analysis of Sox8 and Sox9 during sex determination in the mouse[END_REF][START_REF] Foster | Campomelic dysplasia and autosomal sex reversal caused by mutations in an SRY-related gene[END_REF][START_REF] Qin | Sox9 is sufficient for functional testis development producing fertile male mice in the absence of Sry[END_REF][START_REF] Vidal | Sox9 induces testis development in XX transgenic mice[END_REF][START_REF] Wagner | Autosomal sex reversal and campomelic dysplasia are caused by mutations in and around the SRY-related gene SOX9[END_REF], has a fruit fly ancestor, Sox100B, which was found to be necessary for testis development in Drosophila [START_REF] Nanda | Sox100B, a Drosophila group E Sox-domain gene, is required for somatic testis differentiation[END_REF]. However, the most conserved sex-differentiating factor throughout evolution is DMRT1. Indeed, it has been maintained at the head of the sex determination cascade in reptiles [START_REF] Sun | Dmrt1 is required for primary male sexual differentiation in Chinese soft-shelled turtle Pelodiscus sinensis[END_REF], fishes [START_REF] Matsuda | DMY is a Y-specific DM-domain gene required for male development in the medaka fish[END_REF], and birds [START_REF] Smith | The avian Z-linked gene DMRT1 is required for male sex determination in the chicken[END_REF]. Nevertheless, its functions could have been reduced in mammals since testis differentiates in the absence of DMRT1 (Dmrt1 -/-) in mice [START_REF] Raymond | Dmrt1, a gene related to worm and fly sexual regulators, is required for mammalian testis differentiation[END_REF]. Our results highlight an evolutionary continuum of this gene in testis determination from birds to rabbits and non-rodent mammals in general. Interestingly, even DMRT1 dosage sensibility has been conserved between chicken and rabbits since heterozygous XY DMRT1 +/-male rabbits present secondary infertility with an arrest of spermatogenesis around 2 years of age (data not shown).

DMRT1 position in the rabbit sex-determining cascade

As the early stages of gonadal differentiation in rabbits were not fully characterized, we first determined the expressional profiles of the major sex-determining genes. We observed that DMRT1 expression started at 12 dpc, at the early formation of genital crests, and it was first expressed in the somatic lineage of both sexes, as in mice [START_REF] Lei | Sex-specific differences in mouse DMRT1 expression are both cell type-and stage-dependent during gonad development[END_REF][START_REF] Raymond | Expression of dmrt1 in the genital ridge of mouse and chicken embryos suggests a role in vertebrate sexual development[END_REF] or in humans [START_REF] Garcia-Alonso | Single-cell roadmap of human gonadal development[END_REF]. In the human fetal testis, DMRT1 expression is co-detected with SRY in early supporting gonadal cells, which become Sertoli cells following the activation of SOX9 expression [START_REF] Garcia-Alonso | Single-cell roadmap of human gonadal development[END_REF]. In mice, the Dmrt1 expression starts at E10.5 in both somatic and germinal compartments. However, we showed that germline expression was shifted by 6-8 days compared to the somatic compartment in the rabbit male and female gonads, respectively. These differences are strongly related to the timing of gonadal development in rabbits -which is longer than in mice -and therefore allows better visualization of the different processes. These sequential DMRT1 up-regulations according to cell type and sex also argue in favor of distinct DMRT1 promoters as already described in rats [START_REF] Lei | Distinct transcriptional mechanisms direct expression of the rat Dmrt1 promoter in sertoli cells and germ cells of transgenic mice[END_REF]. For the somatic XY compartment, SOX9 expression appears at 14 dpc in cells expressing both DMRT1 and SRY, suggesting that both factors are required for SOX9 up-regulation. This led to the Sertoli cell differentiation and testicular cords formation from 15 dpc. In the developing ovary, we showed that FOXL2 increases when DMRT1 expression starts to shift from somatic cells to germ cells. Moreover, our results suggested DMRT1 involvement in RSPO1 up-regulation in the ovary.

DMRT1 is required for testis determination in rabbits

In recent years, the advent of new genome editing technologies has made it possible to explore other animal models, such as the goat [START_REF] Boulanger | FOXL2 Is a Female Sex-Determining Gene in the Goat[END_REF] or the rabbit [START_REF] Jolivet | Fetal estrogens are not involved in sex determination but critical for early ovarian differentiation in rabbits[END_REF], and enriching our knowledge on the conservation of ancestral genetic mechanisms in non-rodent mammals. In rabbits, the CRISPR-Cas9 technology allowed us to generate a null mutation of the DMRT1 gene, leading to an absence of detectable protein at homozygosity. Thanks to this model, we could demonstrate that DMRT1 kept its leadership in sex determination also in mammals, where SRY stays the 'switch-on factor' for testis determination, as previously demonstrated in rabbits [START_REF] Song | Mutation of the Sp1 binding site in the 5′ flanking region of SRY causes sex reversal in rabbits[END_REF]. Very early in fetal life, XY fetuses expressing SRY but lacking DMRT1 (DMRT1 -/-) presented a male-to-female sex reversal. Although SRY expression was maintained in XY homozygous mutant gonads, the activation of SOX9 expression was weak in the absence of DMRT1. Accordingly, a few cells expressing SOX9 protein were detectable, but SOX9 target genes expression were not activated in XY DMRT1 -/-gonads. Thus, DMRT1 seems to be required for SRY action on its targets (i.e., SOX9 gene activation) but also for SOX9 functions in the early fetal gonad. Interestingly, a recent study proposed that DMRT1 can act as a SOX9 pioneer factor in the post-natal testis for Sertoli cell identity maintenance [START_REF] Lindeman | The conserved sex regulator DMRT1 recruits SOX9 in sexual cell fate reprogramming[END_REF]. In rabbits, DMRT1 is required for SOX9 and SRY functions, and we hypothesize that DMRT1 might be a pioneer factor for both. In the differentiating genital crest, DMRT1 would be required to increase chromatin accessibility on specific sex-related regions, allowing SRY to bind and activate its targets and particularly the expression of SOX9. The crucial region for SRY binding was identified in mice more than 500 kb upstream of the Sox9 transcription start site and named Enhancer 13 [START_REF] Gonen | Sex reversal following deletion of a single far upstream enhancer of Sox9[END_REF]. Conservation studies identified the homolog of Enhancer 13 in many mammalian species, including humans, cows, and rabbits, and DMRT1 consensus sites were predicted in all mammals examined except mice and rats [START_REF] Gonen | Sex reversal following deletion of a single far upstream enhancer of Sox9[END_REF]. In non-rodent mammals, DMRT1 might be required for chromatin remodeling on the Enhancer 13 region to enable SRY binding and SOX9 expression since the beginning of testis differentiation. In the mouse, which evolved more rapidly, DMRT1 would no longer be necessary for SRY action because the chromatin state of the fetal supporting cells would be more permissive. This could also explain why DMRT1 does not exert any critical function in the fetal testis in mice [START_REF] Raymond | Dmrt1, a gene related to worm and fly sexual regulators, is required for mammalian testis differentiation[END_REF]. In contrast, it is required for the action of SOX9 in the post-natal testis [START_REF] Lindeman | The conserved sex regulator DMRT1 recruits SOX9 in sexual cell fate reprogramming[END_REF], where a sex-specific epigenetic signature was observed [START_REF] Garcia-Moreno | Gonadal supporting cells acquire sex-specific chromatin landscapes during mammalian sex determination[END_REF].

DMRT1 is required for germ cell meiosis and female fertility

In addition to its functions in testis differentiation, DMRT1 also plays a crucial role in the female gonad. Indeed, germ cells did not undergo meiosis in DMRT1 -/-ovaries, and in the absence of oocyte I, germ cell cysts do not break, compromising follicle formation and female fertility. This specific phenotype is highly similar to those observed in ZW chicken ovaries lacking DMRT1 [START_REF] Ioannidis | Primary sex determination in birds depends on DMRT1 dosage, but gonadal sex does not determine adult secondary sex characteristics[END_REF], but is quite different from those described in mice. Even though fewer follicles were observed in Dmrt1 -/-mice ovaries, the female remains fertile [START_REF] Krentz | DMRT1 promotes oogenesis by transcriptional activation of Stra8 in the mammalian fetal ovary[END_REF]. Interestingly in humans, one case involving DMRT1 in premature ovarian failure has been reported [START_REF] Bartels | Normal intelligence and premature ovarian failure in an adult female with a 7.6 Mb de novo terminal deletion of chromosome 9p[END_REF].

In rabbit fetuses, DMRT1 expression was first detected in differentiating ovarian somatic cells, at least until FOXL2 up-regulation. However, DMRT1 has also been observed in fetal germ cells from 20 dpc until meiosis proceeded after birth. Consequently, germ cell pre-meiotic arrest in DMRT1 -/-XX gonads could result from DMRT1 loss-of-function in the germinal or the somatic compartment or both. In the somatic compartment, the absence of DMRT1 in XX homozygous mutants did not seem to disturb the first steps of ovarian differentiation. Nevertheless, deep-sequencing transcriptomics revealed the dysregulated expression of a few genes involved in the WNT/beta-catenin pathway. In particular, RSPO1, a positive regulator of the WNT signaling, was reduced, and DKK1, a negative regulator, was increased (Supplementary file 1 and Supplementary file 2). These two events could have the effect of limiting the beta-catenin action in both somatic and germinal ovarian cells at the beginning of their differentiation. This pathway has proven to be crucial in mice to promote germ cell meiosis (Le [START_REF] Rolle | Arrest of WNT/β-catenin signaling enables the transition from pluripotent to differentiated germ cells in mouse ovaries[END_REF]. Nevertheless, it cannot be the main event explaining the pre-meiotic failure, and the functions of DMRT1 in germ cells are more certainly involved. In mice, DMRT1 was shown to be involved in Stra8 up-regulation in female germ cells and was thus related to the meiotic process [START_REF] Krentz | DMRT1 promotes oogenesis by transcriptional activation of Stra8 in the mammalian fetal ovary[END_REF]. This regulatory action also seems to be done in close collaboration with the retinoic acid pathway [START_REF] Feng | Identification of regulatory elements required for Stra8 expression in fetal ovarian germ cells of the mouse[END_REF]. The sole action of DMRT1 on STRA8 up-regulation cannot explain the phenotype observed in rabbits where the germ cell seems to be unable to leave their pluripotency stage. It has also been demonstrated in male mice that DMRT1 acts as a regulator of spermatogonia pluripotency by directly regulating different pluripotency-associated genes, including POU5F1 [START_REF] Krentz | The DM domain protein DMRT1 is a dose-sensitive regulator of fetal germ cell proliferation and pluripotency[END_REF][START_REF] Zhang | DMRT1 Is Required for Mouse Spermatogonial Stem Cell Maintenance and Replenishment[END_REF]. This path is under exploration in our model in order to try to decipher further the critical role of DMRT1 in the germ line of both sexes.

Materials and methods

Continued on next page

Key resources table

Animals

New Zealand rabbits (NZ1777, Hypharm, Roussay, France) were bred at the SAAJ rabbit facility (Jouyen-Josas, France). All experiments were performed with the approval of the French Ministry MENESR (accreditation number APAFIS#685 and #21451) and following the guidelines issued by the local committee for ethics in animal experimentation (COMETHEA, Jouy-en-Josas). All scientists working directly with the animals possessed an animal experimentation license delivered by the French veterinary services. Hormonal superovulation treatments and surgical embryo transfer procedures were performed as previously described [START_REF] Peyny | Investigating the role of BCAR4 in ovarian physiology and female fertility by genome editing in rabbit[END_REF].

Generation of mutant rabbits

Two guide RNAs were designed (http://crispor.trefor.net/) to target the third exon, as shown in Figure 4-figure supplement 1A. Embryos produced from superovulated females were injected at the single-cell stage with a mixture of the two sgRNAs (10 ng/µl each) and the Cas9mRNA (10 ng/µl) in the injection buffer. Injected embryos were implanted 3-4 hr after into the oviducts of anesthetized recipient rabbits via laparotomy. Details concerning the handling of females and embryos have been described elsewhere [START_REF] Peyny | Investigating the role of BCAR4 in ovarian physiology and female fertility by genome editing in rabbit[END_REF].

Offspring were screened for the presence of InDel mutations using genomic DNA extracted from ear clips [START_REF] Jolivet | Induction of body weight loss through RNAi-knockdown of APOBEC1 gene expression in transgenic rabbits[END_REF]. Founders were detected by PCR using one set of primers (Table 1) surrounding the position of the targeted region in exon III (Figure 4-figure supplement 1A). The amplified fragment was sequenced (Eurofins Genomics, Courtaboeuf, France), and the mutation was deduced by comparing it with the sequence of a wild-type rabbit. The same set of primers was used for the routine screening of descendants. The presence/absence of the Y chromosome was deduced from the amplification of the SRY gene through PCR analyses (Table 1). In the present paper, mentions of the XY or XX genotype always refer to the PCR determination.

XY and XX DMRT1 +/-rabbits were viable until adulthood and did not appear to have any diseases. DMRT1 -/-mutants were obtained by crossing XY DMRT1 +/-and XX DMRT1 +/-animals.

Histological and immunohistological analyses

Immediately after sampling, whole embryos or gonads were immersed in 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS) or Bouin's fixative. After 72-96 hr of fixation at 4°C, tissues were washed three times with PBS, and stored at 4°C in 70% ethanol until paraffin inclusions. Adjacent sections of 5 µm thick were processed using a microtome (Leica RM2245) and organized on Superfrost Plus Slides (J18000AMNZ, Epredia). Before staining or experiments, sections were deparaffinized and rehydrated in successive baths of xylene and ethanol at room temperature.

Hematoxylin-eosin-saffron (HES) staining was performed by the @Bridge platform (INRAE, Jouyen-Josas, France) using an automatic Varistain Slide Stainer (Thermo Fisher Scientific).

ISH was performed using the RNAscope ISH methodology (ACD, Bio-Techne SAS, Rennes, France) when no reliable antibody could be used to characterize the target protein. Briefly, 5 µm sections from PFA-fixed tissue were labeled using RNAscope 2.5HD assay-brown kit (322310, ACD) and 1000 nucleotides long probes designed and produced by the manufacturer (list of all synthesized probes used in Table 2). Brown labeling was observed as a visible signal, and hybridization was considered to be positive when at least one dot was observed in a cell. Immunohistochemistry (IHC) was performed using the ABC amplification signal kit (PK-6100, Vector Laboratories) and DAB enzymatic reaction (SK-4100, Vector Laboratories). Briefly, the antigenic sites were unmasked with a citrate buffer (pH 6; H-3300, Vector Laboratories), and endogenous peroxidases were blocked with a 3% H 2 O 2 solution (H1009, Sigma-Aldrich). Sections were then permeabilized with 1× PBS, 1% bovine serum albumin (A7906, Sigma-Aldrich), and 0.2% saponin (7395, Merck) and incubated overnight at 4°C with primary antibodies (Table 3). Following PBS washes, sections were incubated with biotinylated secondary antibodies (Table 3). After ABC kit incubation and DAB revelation, hematoxylin staining was briefly performed to visualize the whole tissue.

Immunofluorescence (IF) was performed using Tyramide SuperBoost kit for primary rabbit antibody (B40944, Thermo Fisher) as recommended by the manufacturer. Other secondary antibodies used are listed in Table 3.

All stained sections were scanned using a 3DHISTECH panoramic scanner at the @Bridge platform (INRAE, Jouy-en-Josas, France).

Total RNA extraction and RT-qPCR

Immediately after sampling, 16-20 dpc gonads were snap-frozen in liquid nitrogen and stored at -80°C until extraction. Total RNAs were isolated using Trizol reagent (15596018, Life Technologies), Reverse transcription of 50-100 ng RNAs using the Maxima First-Strand cDNA Synthesis Kit (K1641, Thermo Scientific) was down. qPCR with diluted cDNA was performed in duplicate for all tested genes with the Step One system (Applied Biosystems) and Fast SYBR Green Master Mix (4385612, Applied Biosystems). H2AFX and YWHAZ or SF1 (Splicing Factor 1) were used as the reference genes to normalize the results with qBase + software (Biogazelle NV, Ghent, Belgium). The sequences of the primers used are listed in Table 1.

For each experiment, values were plotted using GraphPad Prism Software (GraphPad Software Inc, La Jolla, CA, USA). Statistical analyses of data from 20 dpc control and DMRT1 -/-gonads were performed under R studio software. Because of the small number of samples in each group, comparisons were made using the Kruskal-Wallis rank sum test followed by pairwise permutation t-tests (1000 permutations, p-value adjusted with the Benjamini-Hochberg method).

Nuclear proteins extraction and western blot

Gonads from newborns (1-3 days post-partum) rabbits were collected and snap-frozen in liquid nitrogen and then stored at -80°C. Frozen gonads were crushed in liquid nitrogen using a mortar. Powdered tissue samples were immediately resuspended in homogenization buffer (10 mM HEPES pH 7.7; 25 mM KCl; 2 mM Sucrose; 0.5 mM EGTA pH 8; 0.15 mM Spermin; 0.5 mM Spermidin; 0.5 mM Dithiothreitol (DTT); 2 mM Benzamidin; 0.5 mM Phenylmethylsulfonyl fluoride (PMSF); cOmplete, Mini, EDTA-free Protease Inhibitor Cocktail (Roche, 118361700001)) plus 0.3% IGEPAL (3021, Sigma-Aldrich). After centrifugation for 15 min at 4°C and 3500 rpm, supernatants containing cytosolic proteins were stored at -80°C. Pellets were centrifugated for 1 hr at 4°C and 12,700 rpm and then resuspended with [C-NaCl] buffer (20 mM HEPES pH 7.7; 1.5 mM MgCl 2 ; 0.2 mM EDTA; 25% glycerol; 0.5 mM PMSF; 0.5 mM DTT; 2 mM Benzamidin; cOmplete, Mini, EDTA-free Protease Inhibitor Cocktail). NaCl was added, lysates were rotated for 1 hr at 4°C and then centrifugated for 30 min at 4°C and 12,700 rpm. Supernatants containing nuclear extracts were collected, and the amount of protein was determined by the Bradford method.

Protein (20 µg of each sample) was separated on 4-15% polyacrylamide gel (456-1083, Bio-Rad) and then transferred into a polyvinylidene difluoride membrane. The membrane was blocked in 4% milk (Difco Skim Milk #232100 diluted in PBS-Tween 2%) and incubated overnight with primary antibodies mouse anti-DMRT1 (Table 3) or mouse anti-beta Actin (1/5000; GTX26276, Genetex). After washes, the membrane was incubated for 1 hr at room temperature with the secondary antibody anti-mouse IgG peroxidase-conjugated (1/500; A5906, Sigma-Aldrich). The revelation was performed using Pierce ECL Plus Western Blotting Substrate (32312, Thermo Fisher), and the signal was observed with the Chemi-Doc Touch Imaging System (Bio-Rad). For rehybridization, the membrane was stripped for 10 min in Restore Western Blot stripping buffer (21059, Thermo Fisher).

RNA-sequencing and bioinformatics analysis

Total RNAs were extracted from control and DMRT1 -/-rabbit gonads at 20 dpc (n = 3 for each phenotype and each sex). Total RNA quality was verified on an Agilent 2100 Bioanalyser (Matriks, Norway), and samples with a RIN >9 were made available for RNA-sequencing. This work benefited from the facilities and expertise of the I2BC High-throughput Sequencing Platform (https://www.i2bc.parissaclay.fr/sequencing/ng-sequencing/ Université Paris-Saclay, Gif-sur-Yvette, France) for oriented library preparation (Illumina Truseq RNA Sample Preparation Kit) and sequencing (paired-end 50-35 bp; NextSeq500). More than 37 million 50-35 bp paired-end reads per sample were generated. Demultiplexing was done (bcl2fastq2-2.18.12), and adapters were removed (Cutadapt1.15) at the I2BC High-throughput Sequencing Platform. Only reads longer than 10 pb were used for analysis. Quality control of raw RNA-Seq data was processed by FastQC v0.11.5. Reads were mapped on all the genes of a better-annotated rabbit genome. Indeed, we improved the current reference rabbit transcriptome (OryCun2.0; Oryctolagus cuniculus, Ensembl version 106). For this purpose, we extended the 5′ and 3′-UTRs of genes using rabbit gonad RNA-seq data available in public databases (https://www.ncbi.nlm.nih.gov/bioproject/PRJEB26840). In addition, the annotation and for some of them, sequences of 22 marker genes of gonadal differentiation missing or wrong in OryCun2.0 was added or fixed to this genome assembly. Then, after mapping with STAR version 2.5.1b [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF], reads were counted using FeatureCounts version 1.4.5 [START_REF] Liao | featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[END_REF]. Data normalization and single-gene level analyses of differential expression were performed using DESeq2 [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]. Differences were considered to be significant for Benjamini-Hochberg adjusted p-values <0.05, and absolute fold log2FC >1 [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. RNA-seq data were deposited via the SRA Submission portal (https://www.ncbi.nlm.nih.gov/sra/PRJNA899447), BioProject ID PRJNA899447.

Figure 1 .

 1 Figure 1. SRY, DMRT1, and POU5F1 location during early gonadal development. (A) Key stages of gonadal development in rabbits with 31 days of gestation. Germ cells are first detected at 9 days post-coïtum (dpc), before the genital ridge formation, which occurs between 10 and 12 dpc. In XY gonads, testicular cords begin forming at 16 dpc, and germ cells enter meiosis a few months after birth. In XX gonads, the ovigerous cords appear at 20 dpc, and meiosis begins around birth. Location of SRY, DMRT1, and POU5F1 by in situ hybridization (RNAscope technology) on XY and XX control gonads at (B) 12 dpc or (C) 13 dpc. Dotted line: developing genital crests. Yellow arrowheads: coelomic epithelial cells expressing DMRT1. Scale bar = 50 µm.

Figure 2 .

 2 Figure 2. Somatic markers location during testis differentiation. Location of SRY by in situ hybridization (RNAscope technology), DMRT1, and SOX9 by immunohistochemistry on XY control testes from 14 to 18 dpc. The dotted line at 15 dpc: territory with cells expressing SRY and DMRT1 but not SOX9. Yellow arrowheads: tunica albuginea in formation. Scale bar = 50 µm. The online version of this article includes the following figure supplement(s) for figure 2: Figure supplement 1. DMRT1 and SOX9 co-location on 14 dpc XY control gonad.

Figure supplement 2 .

 2 Figure supplement 2. Identification of PAX8-positive cells in 15 dpc control gonads.

Figure 3 .

 3 Figure 3. Somatic markers location and expression during ovarian differentiation. (A) Immunostaining of DMRT1 on XX control ovaries from 14 to 18 dpc. (B) Quantitative RT-PCR (RT-qPCR) analyses of RSPO1, DMRT1, and FOXL2 expression from 16 to 20 dpc in control gonads of both sexes. The error bars correspond to the standard error of the mean (n=3-5) (C) RSPO1 in situ hybridization (RNAscope technology), immunostaining of DMRT1 and FOXL2 on 20 dpc control ovaries. Scale bar = 50 µm. The online version of this article includes the following figure supplement(s) for figure 3:

Figure supplement 1 .

 1 Figure supplement 1. DMRT1 and POU5F1 co-detection in control gonads.

Figure 4 .

 4 Figure 4. Ovarian-like morphology and transcriptomic signature of XY DMRT1 -/-gonads at 20 dpc. (A) Hematoxylin and eosin staining of gonads sections from control and DMRT1 -/-20 dpc rabbits. The enlarged area shows the characteristic ovarian surface epithelium found on XY DMRT1 -/- gonads. Scale bar = 50 µm. Heatmap representation of (B) 3460 deregulated genes (adjusted p-value <0.05 and |log2FC| > 1) or (C) 27 selected genes between XY control, XY DMRT1 -/-, XX DMRT1 -/-, and XX control at 20 dpc. The online version of this article includes the following figure supplement(s) for figure 4: Figure supplement 1. DMRT1 mutation using CRISPR/Cas9 in rabbits.

DujardinFigure 5 .

 5 Figure 5. Somatic markers expression and location on control and DMRT1 -/-gonads at 20 dpc. Quantitative RT-PCR (RT-qPCR) analyses of (A) testicularrelated differentiation genes (SOX9, AMH, DHH, and SRY) or (B) ovarian-related differentiation genes (FOXL2, CYP19A1, and RSPO1) in XY control, XY DMRT1 -/-, XX DMRT1 -/-, and XX control gonads (n = 4-5) at 20 dpc. Statistical analyses were performed using the non-parametric Kruskal-Wallis test, followed by a pairwise permutation test: *p-value <0.05; ns: non-significant. (C) Immunostaining of DMRT1, SOX9, and FOXL2 on XY control, XY DMRT1 -/-, XX DMRT1 -/-, and XX control gonad sections at 20 dpc. Scale bar = 50 µm.

Figure 6 .

 6 Figure 6. Evolution of gonadal morphogenesis in XY and XX DMRT1 -/-rabbits. Hematoxylin and eosin staining of gonad sections from XY and XX DMRT1 -/-gonads and XX control ovaries at 3 days post-partum (dpp), 18 dpp, and in adulthood (4-9 months). The enlargements for the first two panels correspond to the nuclei pointed by an arrow. PL: preleptotene stage; L: leptotene stage; Z: zygotene stage; D: diplotene stage; F: ovarian follicle; CL: luteal cells. Scale bar = 50 µm. The online version of this article includes the following figure supplement(s) for figure 6:

Figure supplement 1 .

 1 Figure supplement 1. POU5F1 and Ki67 location on control and DMRT1 -/-gonads at 3 dpp.

Figure supplement 2 .

 2 Figure supplement 2. Evolution of gonadal size in XY and XX DMRT1 -/-rabbits.

  Dujardin et al. eLife 2023;12:RP89284. DOI: https://doi.org/10.7554/eLife.89284 14 of 20

Table 1 .

 1 Primers used for genotyping PCR or quantitative RT-PCR (RT-qPCR) analyses.

	Gene	Forward (5′-3′)	Reverse (3′-5′)
	Genotyping PCR
	DMRT1	TTTG AGCT GTGT CCCC AGAG T	ACCT CCCC AGAA GAAG AATC G
	SRY	GTTC GGAG CACT GTAC AGCG	GCGT TCAT GGGT CGCT TGAC
	RT-qPCR analyses
	SRY	TGCT TACA CACC AGCC AAAC A	TTCC TGGC CGCT CACT TTAC
	DMRT1	GGAG CCTC CCAG CACC TTA	TGCA TCCT GTAC TGCG AACT CA
	SOX9	GGCT CCGA CACC GAGA ATAC AC	GAAC TTGT CCTC TTCG CTCT CCTT
	CYP19A1	GGAA GAAT GCAT CGAC TTGA GTT	GGGC CCAA AACC AAAT GGT
	ESR1	TCCT CATC CTCT CCCA CATC	AGCA TCTC CAGC AACA GGTC
	RSPO1	GCCC GCCT GGAT ACTT CGA	GGTG CAGA AGTT GTGG CTGA A
	FOXL2	TTTC CCCT TTCC CCCA TCTG	CTGA ACCT TGCA CCCA GCAT
	AMH	GCTC ATCC CCGA GACC TAC	CATC TTCA ACAG CAGC ACC
	DHH	GCAA TAAG TACG GGCT GCTG	GGCC AGGG AGTT ATCA GCTT
	H2AFX	ACCT GACG GCCG AGAT CCT	CGCC CAGC AGCT TGTT GAG
	YWHAZ	GGGT CTGG CCCT TAAC TTCT CT	AGCA ATGG CTTC ATCA AAAG C
	SF1 (splicing factor 1)	GCTT CCGA CTGC AAAT TCCA	TCAC CCAG TTCA GCCA TGAG

Table 2 .

 2 Synthesized probes used for in situ hybridization.Dujardin

	Gene name	RNAscope probe catalog number	Transcript accession number
	DMRT1	410481	XM_002708188.1
	SRY	803191	AY785433.1
	POU5F1	513271	NM_001099957.1
	RSPO1	488231	XM_002720657.2

et al. eLife 2023;12:RP89284. DOI: https://doi.org/10.7554/eLife.89284 15 of 20 purified with the RNeasy Micro kit (74004, QIAGEN) following the manufacturer's instructions, and then DNAse treated (1023460, QIAGEN). RNAs were quantified with a Qubit Fluorometric Quantification kit (Q32852, Life Technologies).

Table 3 .

 3 List of antibodies used for immunohistochemistry (IHC), immunofluorescence (IF), or western blot (WB).

		Antibody name	Reference	Dilution
				1/500 (IHC); 1/200 (IF); 1/100
		Mouse monoclonal to DMRT1	sc-377167 (Santa Cruz)	(WB)
		Rabbit polyclonal to SOX9	Francis Poulat	1/500 (IHC); 1/200 (IF)
		Rabbit polyclonal to FOXL2	Boulanger et al., 2014	1/500 (IHC)
		Goat polyclonal to POU5F1	sc-8628 (Santa Cruz)	1/500 (IHC); 1/200 (IF)
		Rabbit polyclonal to PAX8	10336-1-AP (Proteintech)	1/6000 (IHC)
	Primary antibodies	Rabbit monoclonal to Ki67	MA5-14520 (Thermo Scientific)	1/500 (IHC)
		Horse anti-rabbit	BA-1100 (Vector Laboratories)	1/200
			Included in M.O.M. kit (BMK-2202, Vector
	Biotinylated secondary	Anti-mouse	Laboratories)	1/200
	antibodies (IHC)	Horse anti-goat	BA-9500 (Vector Laboratories)	1/200
		Poly HRP-conjugated goat anti-rabbit	B40944 (Invitrogen)	No diluted
		DyLight 488 goat anti-mouse	072-03-18-06 (KPL)	1/200
	Secondary antibodies (IF)	Alexa Fluor 594 chicken anti-goat	A21468 (Life technologies)	1/200

Dujardin et al. eLife 2023;12:RP89284. DOI: https://doi.org/10.7554/eLife.89284 16 of 20

  Dujardin et al. eLife 2023;12:RP89284. DOI: https://doi.org/10.7554/eLife.89284 18 of 20The following dataset was generated:

	Author(s)	Year	Dataset title	Dataset URL	Database and Identifier
	Dujardin E, Andre M,	2023	RNA-seq of Oryctolagus	https://www. ncbi.	NCBI Sequence Read
	Dewaele A, Mandon-		cuniculus	nlm. nih. gov/ sra/	Archive, PRJNA899447
	Pépin B, Poulat			PRJNA899447
	F, Frambourg A,
	Thépot D, Jouneau L,
	Jolivet G, Pailhoux E,
	Pannetier M

Acknowledgements

The authors would like to thank Patrice Congar, Gwendoline Morin, and all the staff of the facility (SAAJ, INRAE, Jouy-en-Josas, France) for the care of the rabbits, Erwana Harscoët and Nathalie Daniel for injecting the embryos, Laurent Boulanger for the sgRNAs design, Nathalie Daniel-Carlier for the mutation characterization, Julie Rivière and Marthe Vilotte (UMR GABI, INRAE, Jouy-en-Josas, France) for their assistance on the histological platform (@Bridge platform) and the access to the virtual slide scanner, Namya Mellouk for her contribution to the dissection of the gonads, Simon Herman (Université Paris-Saclay, France) for his contribution to the improvement of the rabbit genome annotation during his master's internship in our team (DGP, UMR BREED, INRAE, Jouy-en-Josas, France), and Andrew Crawford (Academic Writing Center, Centralesupélec, France) for the English proofreading. Francis Poulat kindly provided the SOX9 antibody. We acknowledge the sequencing and bioinformatics expertise of the I2BC High-throughput sequencing facility, supported by France Génomique (funded by the French National Program 'Investissement d'Avenir' ANR-10-INBS-09). We are grateful to the genotoul bioinformatics platform Toulouse Occitanie (Bioinfo Genotoul, France, https://doi. org/10.15454/1.5572369328961167E12) for providing computing and storage resources. This work was supported by ANR (Agence Nationale de la Recherche) grants (RNA-SEX: ANR-19-CE14-0012; ARDIGERM: ANR-20-CE14-0022). ED was supported by the ANR RNA-SEX and the PHASE department of INRAE.

Funding

Funder

Grant reference number Author

Data availability RNA-seq data were deposited via the SRA Submission portal, BioProject ID PRJNA899447.

Additional information

Ethics All experiments were performed with the approval of the French Ministry MENESR (accreditation number APAFIS#685 and #21451) and following the guidelines issued by the local committee for ethics in animal experimentation (COMETHEA, Jouy-en-Josas). All scientists working directly with the animals possessed an animal experimentation license delivered by the French veterinary services.

Peer review material

Additional files

Supplementary files

• Supplementary file 1. List of DEGs in DMRT1 -/-gonads (XY and XX) compared to control gonads (XY and XX) (adjusted p-value <0.05 and |log2FC| > 1). List of deregulated genes (DEGs) between KO-XY vs Control-XY gonads (sheet 1), KO-XY vs KO-XY gonads (sheet 2), KO-XX vs Control-XX gonads (sheet 3), and Control-XY vs Control-XX gonads (sheet 4). The gene name of DEGs was based on their annotation or human homology [START_REF] Craig | Augmented annotation and orthologue analysis for Oryctolagus cuniculus: Better Bunny[END_REF].

• Supplementary file 2. Clustering and expression values (TPM, transcripts per million) of the 3460 DEGs. Cluster membership of the 3460 DEGs with their expression data (TPM) according to the four genotypes (XY control, XY DMRT1 -/-, XX DMRT1 -/-, XX control).

• MDAR checklist


Igor Maingonnat

Gilles Tissot

Noé Lahaye

Correlated structures in a balanced motion interacting with an internal wave

Keywords: Internal tide interactions, spectral proper orthogonal decomposition, broadband proper orthogonal decomposition

Characterizing the loss of coherence of an internal tide propagating through mesoscale turbulence has been a major challenge in oceanography, particularly due to its implications for the interpretation of satellite data. In this paper, we intend to study the correlations between a balanced motion and the incoherent part of a wave in an idealised configuration. We introduce a new modal decomposition technique, named broad-band proper orthogonal decomposition (BBPOD), which consists in performing a proper orthogonal decomposition (POD) on complex demodulated variables. After connecting BBPOD to the standard SPOD, we show that BBPOD, coupled with the extended POD technique enables us to associate the principal components of the incoherent field to the slow flow structures responsible of this loss of coherence through triadic interactions with the incident wave.

Introduction

Internal tides, generated by interactions between the barotropic tide and topographic features such as ridges or continental slopes, are ubiquitous in the ocean, playing a crucial role in vertical mixing and energy transport. They propagate over large distances, encountering regions with energetic mesoscale turbulence, and they lose their fixed phase relationship with the astronomical forcing, a phenomenon known as incoherence. This loss of coherence, highly unpredictable, complicates for example our ability to disentangle internal tides and low-frequency turbulent signals from satellite data [START_REF] Richman | Inferring dynamics from the wavenumber spectra of an eddying global ocean model with embedded tides[END_REF]).

To face these difficulties, there have been various studies to better understand and predict the impact of a balanced (turbulent) jet on the inertia-gravity wave field properties. For instance, [START_REF] Savva | Scattering of internal tides by barotropic quasigeostrophic flows[END_REF] and [START_REF] Ward | Scattering of gravity waves by potential vorticity in a shallow-water fluid[END_REF] studied the internal tide scattering. [START_REF] Ponte | Incoherent signature of internal tides on sea level in idealized numerical simulations[END_REF] examined the incoherence time scales for different turbulent configura-tions, and [START_REF] Dunphy | Low-mode internal tide propagation in a turbulent eddy field[END_REF] quantified the interaction terms via vertical-mode projection of the linearized Boussinesq equations.

We consider in the present paper a data-driven approach to study from idealised numerical simulations, the structures of the jet, which are correlated with the incoherent contributions of the wave field. Extracting coherent structures from a data set can be performed in the spectral domain by spectral proper orthogonal decomposition (SPOD) [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. Some attempts to connect these reduced features with non-linear interactions have been proposed for instance in [START_REF] Karban | Self-similar mechanisms in wall turbulence studied using resolvent analysis[END_REF] by identifying the contribution of the non-linear term correlated with the SPOD mode through extended proper orthogonal decomposition (EPOD, Boree (2003)). Unfortunately, non-linear forcing is a quantity more difficult to interpret and associate with physical mechanisms than a pressure or velocity field for instance. We propose instead a new broad-band proper orthogonal decomposition (BBPOD), derived from the SPOD. By taking advantage of the strong time scale separation between the two dynamics and considering small-amplitude wave, this formulation allows us to connect the non-linear interactions between the slow motions of the jet and the incoherent contributions of the wave through the EPOD method. The non-linear interactions are here understood as triadic interactions with the slow motion that will lead to a broadband frequency structure of the incoherent wave field.

The plan is as follows. We will begin by describing the model used for our simulations in section 2. In section 3.1, the original SPOD method is reviewed and a connection to the proposed BBPOD method is made (section 3.2). Sections 4 and 5 will summarize the study and bring some perspectives.

Model

The propagation of internal tides through a nearly-balanced jet is examined in a one layer rotating shallow water (RSW) model. The equations are non-dimensionalized as follows. The characteristic time-scale is the inertial time T = f -1 , inverse of the Coriolis frequency. The reference length scale L is chosen to be of the order of the jet thickness, taken equal to the first Rossby radius of deformation R d , such that the Burger number B u = R 2 d /L 2 is equal to one. A beta-plane approximation, with parameter β, accounts for the effect of rotation, and a radiative damping term α is added on the continuity equation, (following e.g. [START_REF] Brunet | Empirical Normal Modes versus Empirical Orthogonal Functions for Statistical Prediction[END_REF].

Adequat artificial hyperviscosity, consistent with energy dissipation and conservation of angular momentum, is also used (see [START_REF] Ochoa | Lateral Friction in Reduced-Gravity Models: Parameterizations Consistent with Energy Dissipation and Conservation of Angular Momentum[END_REF].

Let Ω be a bounded subset of R 2 , the equations without viscosity, defined on Ω × R + , are

∂ t h + B u div v + R o ((v•∇)h + h div v) = -αh (1)
∂ t v + (1 + βy)v ⊥ + R o (v•∇)v = -∇h, (2)
where v = (u, v) is the horizontal velocity, v ⊥ = (-v, u) T , h is the sea surface height (SSH)

anomaly and R o = U /(f L) is the Rossby number. -20, 20] discretized on a 128 × 512 grid. The simulation is initialized with an eastward zonal jet at geostrophic equilibrium with a small perturbation to trigger the instabilities. During the whole experiment, a northward propagating plane wave with frequency fixed at ω = 2 is generated in a nudging layer in the South of the domain, and an eastward wind forcing is applied to maintain the balanced current. A spin-up phase allows the jet to evolve toward a statistically stationary state. The experiment continues with snapshots stored every 0.1 wave period, such that a sufficiently long series q = (u, v, h) T of 12000 snapshots are collected, representing nearly 4000 wave periods. Note that a sufficient sampling in time is required for extracting the wave field by filtering in 3.2. An example of snapshot associated with one run of the simulation is shown in Figure 1.

This simple model can sustain waves interacting with a turbulent flow. Northward rotating shallow water waves are inertia gravity waves satisfying the dispersion relation

k 2 y = ω 2 -(1 + βy) 2 B u , (3)
with k y denoting the spatial wavenumber in the meridional direction, which is function of the tidal frequency ω.

Methods

The methods to extract correlated structures in a flow field are described in this section. In the following we denote our state vector by q = (u, v, h) T .

Spectral proper orthogonal decomposition

Spectral proper orthogonal decomposition is an extension of POD methods (see [START_REF] Lumley | The Structure of Inhomogeneous Turbulent Flows[END_REF], and aims at extracting coherent structures in spectral space from numerical or experimental data. The data are assumed to be issued from a statistically stationary random process, which is verified thanks to the hypothesis of small-amplitude waves, such that the balanced flow is marginally impacted by the wave propagation. The modes obtained from this decomposition are space-time uncorrelated from each other.

We interpret our state vector as a zero-mean second order random process q (which can be obtained by subtracting the time averaged field beforehand) indexed over Ω × R + , a subset of spatio-temporal variables. Consider q(ω) the Fourier transform of the process at angular frequency ω. We assume that each realization of q(ω) belongs to

L 2 W (Ω, C 3 ) = {g : Ω → C 3 , Ω g * (z)W(z)g(z)
dz < ∞}, the ensemble of square integrable functions in space (periodic at boundary) relatively to a positive definite weight matrix W: Ω → M n (R). The superscipt • * denotes the transpose-conjugate operation. The matrix W is chosen such that the L 2 -norm approximates the energy of the RSW model ( 1)-(2):

∥q∥ 2 W = 1 2 Ω (1 + λ h)(u 2 + v 2 ) dxdy + 1 2B u Ω h 2 dxdy, (4)
where λ refers to the deviation of the isopycnal taken equal to R o /B u and h corresponds to the temporal-mean of SSH anomaly. The term 1 + λ h is assumed to remain strictly positive (which is necessary for ensuring positiveness of the norm), since we do not consider zero or negative sea-surface height.

Let us recall briefly the basic principle of SPOD. By stationarity the auto-correlation function

C(x, y, t, x ′ , y ′ , t ′ ) = E(q(x, y, t) ⊗ q * (x ′ , y ′ , t ′ )) satisfies : C(x, y, t, x, y, t ′ ) = C(x, y, x ′ , y ′ , t -t ′ ), (5)
where ⊗ denotes the dyadic product of two vectors in R 3 . The objective of SPOD method is therefore to find a deterministic function ψ satisfying the Fredholm equation :

Ω S(x, y, x ′ , y ′ , ω)W(x ′ , y ′ )ψ(x ′ , y ′ , ω) dx ′ dy ′ = λ(ω)ψ(x, y, ω), (6)
where S(x, y, x ′ , y ′ , ω) = R C(x, y, x ′ , y ′ , τ )e -iωτ dτ is the cross spectral density matrix (CSD). The Karhunen-Loeve theroem (Loève 1955) states that ∀ω ∈ R, equation ( 6) has an infinite number of solutions (λ j (ω), ψ j (ω)) ∞ j=1 such that: (ψ j (ω)) ∞ j=1 is an orthonormal basis in L 2 W (Ω, C 3 ) where we can expand the Fourier transform of the field into structures uncorrelated from each other:

q(x, y, ω) = ∞ j=1 a j (ω)ψ j (x, y, ω) (7) with E(a j (ω)a j ′ (ω ′ )) = λ j (ω)δ j,j ′ (ω -ω ′
) and ∀j, λ j (ω) ≥ 0. Moreover, a truncated expansion at order n will maximize the mean energy (defined by the norm ∥ • ∥ 2 W ), compared to any other decomposition of the same order, and SPOD provides an optimal decomposition of the CSD:

E(∥ q(ω)∥ 2 W ) = ∞ j=1 λ j (ω) (8) S(x, y, x ′ , y ′ , ω) = ∞ j=1 λ j (ω)ψ j (x, y, ω)ψ * j (x ′ , y ′ , ω). (9)
Therefore, by solving the Fredholm equation at the tidal frequency, we are able to separate the fast from the slow component of the process by a spectral decomposition, expanded onto an orthonormal basis optimal energetically.

It can be remarked that when there is an homogeneous direction, it is possible to compute the Fredholm equation in spectral space with respect to this direction. Each wave-number can be computed independently, and each corresponding Fourier component is solution of the original Fredholm problem. For more complete description of SPOD (see [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF][START_REF] Schmidt | Guide to Spectral Proper Orthogonal Decomposition[END_REF]. In our case, the domain is homogeneous in the x direction, but we still consider the 2D problem for the physical relevance of the extended BBPOD problem presented in section (3.2).

Broadband proper orthogonal decomposition

The BBPOD consists in estimating by complex-demodulation (see [START_REF] Godfrey | An exploratory study of the bi-spectrum of economic time series[END_REF] the correlation tensor of a band-passed-filtered signal, and its eigenfunctions. The associated algorithm is presented in this section. It will be compared with the Welch's method used in SPOD (see [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF][START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]) and it will be shown that under some hypotheses both algorithms are equivalent. It can be noted that in [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF], the connection with complex demodulation is briefly mentioned, and it is leveraged here with the computation of the eigenvectors and eigenvalues of the correlation tensor for a POD decomposition. We highlight by this equivalence, that due to windowing in the Welch's method, the features extracted by SPOD possess a spectral component with a thick frequency band. For the Broadband POD algorithm, this frequency band is explicit, and chosen through the definition of a filter. In the RSW model, the incoherent wave field possesses a broad band structure and Broadband POD allows us to obtain a complete decomposition of this field.

Complex demodulation of the wave field

Given a temporal series x t , the principle of complex demodulation consists in the following computation :

x d = L(x t e -iωt ), (10)
where x d is called the complex demodulated signal, L is a low-pass filter and ω is the demodulated frequency. Complex demodulation enables to extract the slowly varying amplitude envelop and phase deviations of the wave field. These variations are associated with wave incoherence.

First, we decompose the process q(x, y, t) in a sum of q j (x, y, t) and q ω (x, y, t), representing the jet and wave contributions, respectively. The wave field can be expressed as:

q ω (x, y, t) = q d (x, y, t)e iωt , ( 11
)
where ω is the tidal frequency and q d is a complex field slowly varying in time, accounting for the incoherence. The latter is then extracted by the complex demodulation of q (10), with a filter designed such that it isolates the slow variations of the background flow. As a consequence, through the e -iωt shift, the spectral-band of the extracted wave contains all effects of triadic interactions between the slow motions and the coherent wave. The jet contribution q j is obtained by directly filtering q (zero frequency), including here the time-average. In this study, we assume that q is dominated by the superimposition of q j and q ω , which interact non-linearly. Note that for simplicity, a uniform time sampling of the data was chosen but this can be adapted to a non-uniform one, e.g. by interpolation on a regular grid or using filtering algorithms to irregular sampling.

Link with SPOD

For making the link with SPOD to estimate the CSD of a signal x t , we write the temporal filter in a discrete form, with a time spacing ∆t and t k = k∆t. A wide class of linear filter can be expressed as the convolution:

L(x t ) = m i=-m b i x t-i , (12)
where (b i ) -m≤i≤m are the discrete coefficients of the filter. Then,

(L(x t e -iωt )) j = m k=-m b k x j-k e -iωt j-k (13) = j+m k=j-m b j-k x k e -iωt k . ( 14
)
In the Welch method, x t is subdivised into possibly overlapping blocks of size N and with an overlap N o . A Fast Fourier Transform is performed on each windowed block to extract the Fourier component at the tidal frequency, denoted X l ω where l is the block index. We define

X l ω = N k=0 x k+l(N -No) W k e -iωt k , (15)
where W k is a window function. By changing variable k ′ = k + l(N -N o ), we obtain

X l ω = e iωt l(N -No ) N +l(N -No) k ′ =l(N -No) x k ′ W k ′ -l(N -No) e -iωt k ′ . ( 16
)
Assuming that the window function is symmetrical in the middle of each block (which is verified for most windows used in the literature), i.e W k = W N -k , equation ( 16) gives

X l ω = e iωt l(N -No) N +l(N -No) k ′ =l(N -No) W N +l(N -No)-k ′ x k ′ e -iωt k ′ . (17)
Finally, by choosing W k = b k which sets m = N 2 , relation ( 14) yields

X l ω = e iωt l(N -No ) (L(x t e -iωt )) N 2 +l(N -No) . (18)
Consequently, up to a phase, the Welch method corresponds to the computation of the complex demodulated signal sampled every N -N o . The phase shift cancels when computing the correlation over N b blocks :

N b l=0 X l ω X l ω * = N b l=0 (L(x t e -iωt )) N 2 +l(N -No) (L * (x t e -iωt )) N 2 +l(N -No) . (19)
The CSD can thus be obtained by computing the correlation over N b snapshots sampled every N -N o of the complex demodulated signal, which is proving the equivalence between the two methods for appropriate numerical parameters.

As said before, it can be remarked that the window function acts as a filter in the Welch procedure, but without giving an explicit expression of the frequency band. Moreover, if we aim at studying the whole broadband field, the classical SPOD (with its single frequency interpretation when the window function is designed to only estimate a pure harmonic) requires the computation of a set of spatial modes at each discrete frequency in the peak. In comparison, BBPOD computes at only one frequency the dominant modes of the wave field. This has the drawback for the SPOD algorithm to include at each frequency a few modes of nearby frequency due to the convolution with a finite length window, leading to misinterpretation or counting the same mode several times in the reconstruction of the whole wave field. More precisely, BBPOD gathers the first SPOD mode of each frequency providing a simpler representation of the whole wave field and avoids spurious modes in the reconstruction. Another important remark is that the connection of our method with SPOD lies on the fact that the spectral peak is sufficiently narrow corresponding of a sufficiently large window in the Welch method guaranteeing a good estimation of the CSD.

The final step is to compute the eigenvalues and eigenvectors of the CSD, estimated by complex demodulation. The complex demodulated signal at frequency ω is then decomposed on this BBPOD basis:

q d (x, y, t) = j a j (ω, t)ψ j (x, y, ω), (20)
with a j (ω, t) = Ω q * d (x, y, t)W(x, y)ψ j (x, y, ω) dxdy slowly varying coefficients.

Extended broadband proper orthogonal decomposition

Non-linear interactions between slow variations of the background flow and the wave induce incoherent wave contributions through triadic interactions. A major interest of BBPOD is that it allows us to study the correlations between the slow perturbations of the jet, and the incoherence of the fast wave field, by extracting the slowly varying complex demodulated amplitudes. To that end, we propose to apply the concept of EPOD in the framework of BBPOD. EPOD, originally presented in [START_REF] Boree | Extended proper orthogonal decomposition: A tool to analyse correlated events in turbulent flows[END_REF], enables to identify the part of a target field correlated with a given POD mode. The target field we first consider is the balanced motion obtained by low-pass filtering. The p-th extended POD mode of the slow motion q j correlated to the p-th BBPOD coefficient of the wave a ω p is defined by :

χ j p = E(q j a p (ω)) λ p (ω) , (21)
where the expectation operator is a temporal average over snapshots. The mode χ j p will be referred to as direct EPOD, and a p (ω)χ j p represents the part of the jet correlated with the p-th broadband POD coefficient of the wave.

Complementary to direct EPOD, we define as well an inverse EPOD, applying BBPOD for ω = 0, thus obtaining an orthonormal basis representative of the jet variability, and considering the complex demodulated of the wave contribution q d as the target field. In this situation, we identify the contribution of the incoherent wave field correlated with the BBPOD coefficients of the jet noted a p (ω = 0). Therefore, in the following direct EPOD refer to jet modes while inverse EPOD refer to wave modes.

Results

This section details the numerical results carried out in this work. The goal is to understand the non-linear interactions between the balanced motion and the wave field by means of the methods presented in section 3, that will extract the correlations between the two dynamics.

Figure 2 shows the real part of the first three energy-scaled BBPOD modes at the tidal frequency, sorted by decreasing energy (4). Specifically here, we plot the scaled quantities λ j ψ j to highlight their respective energy contribution, but the analysis was performed with the normalized modes. To define the low pass filter L, a fourth order Butterworth filter has been taken with a frequency cut-off equals to the typical frequencies of the jet, such that the whole wave-field scattered by the jet is captured. The first mode containing the most energy corresponds to the coherent wave field, that is defined here as the part of the wave correlated to the tidal forcing and thus phase-locked in time, since it corresponds to the time-average of the complex-demodulated signal. Other modes account then for the incoherent structures as space-time decorrelated of the mean. They show deflections in opposite direction, with mode-number m x = ±2, and with meridional mode-number m y expected to be determined by the dispersion relation (3) and the structure of the jet (see Bühler 2014)). This is a consequence of the homogeneity in the zonal direction, guaranteed by our idealised set-up, which makes SPOD modes equal to Fourier modes in x.

Table 1 shows the normalised energy contained in the first five BBPOD modes and in the incoherence field. It is computed for the domain without sponge regions noted Ω = [0, 20] × [-12, 12] and for the subset Ω N = [0, 20] × [-2, 12] representing the region of incoherences.

In the total domain, the coherent plane wave accounts for 57.35 % of the energy and the first three modes 81 % of the energy. This shows that linear effects dominate the wave propagation, because we are in a configuration for which incoherence is relatively weak. As explained in [START_REF] Ward | Scattering of gravity waves by potential vorticity in a shallow-water fluid[END_REF], a more energetic jet interacting with a tide with higher frequency would result in a more energetic incoherent wave field for a RSW simulation. As expected, incoherent modes have nearly all their energy in the upper part while it is rather equally distributed for the coherent mode. In this region, incoherences are dominant due to non-linearities in the center of the domain, which leads to an increase of the amplitude of the modes with y, even though the coherent mode still contains the most energy. Figure 3 shows the energy contained in the reconstructed incoherent wave field in terms of cumulative energy. The first two incoherent modes account for 58 % of the total incoherent energy. For 6 modes, 80 % of the energy is A SPOD decomposition was also performed in [START_REF] Egbert | An Approach to Empirical Mapping of Incoherent Internal Tides With Altimetry Data[END_REF] for a realistic HY-COM simulation. Similar results were found for the energy captured by the SPOD modes at the M2 frequency, which is similar to our forcing frequency. This encourages to extend this analysis to a more realistic simulation. showing that the part of the jet correlated with the coherent field is the stationary mean flow, which connects the first order statistic of both dynamics. This results is the expression of the linear propagation of the incident wave through the mean flow (in time and in the x-direction).

The second and third EPOD modes represent the meanders of the jet, consisting of a vortex train. The modes are defined up to a complex constant. We can see that χ j 2 is approximately proportional to χ j 3 , suggesting that meandering of the jet generates eastward and westward propagating incoherent perturbations.

The opposite procedure is done next by inverse EPOD. Figure 5 shows the leading three jet (energy-scaled) BBPOD modes associated with ω = 0 and the associated EPOD modes of the wave. As in direct EPOD, the first mode is showing the correlation between the mean of both dynamics. For the second and third modes, BBPOD shows jet meandering structures identical to the direct EPOD modes (up to a phase). However, their associated inverse EPOD indicate a standing wave corresponding to the sum of the two direct broadband POD modes that represent left and right deflections, expressing the fact that the EPOD and BBPOD association is not bijective. A given meander will give rise to a wave deflected eastward nearly as much as a wave deflected westward, and their relative contribution is not possible to predict, but the present result suggests that an approximately equal repartition would constitute a fair estimate.

Globally, we can infer that the most energetic triadic interactions between the jet and the plane wave are caused by meandering structures of the jet generating an incoherent field of the form of a standing wave in the zonal direction x and propagating northward. The rising of a standing wave suggests in average an equal distribution of the eastward and westward propagating modes. This gives nice perspectives of incorporating this methodology in the context of data assimilation to estimate the incoherent wave field with some knowledge of the slow balanced flow through triadic interactions. λ j ψ j .

Summary and Perspectives

This study investigates the correlated processes at stake within the non-linear interactions of a plane wave propagating through a meso-scale oceanic current. For that purpose, the method of broadband POD is introduced and connected to extended POD. A first methodological point is done to show that BBPOD is equivalent to the common Welch method to estimate the cross spectral density matrix, to which SPOD are the eigenfunctions. More importantly for our study, connecting these two methods allows us to understand that the most energetic features of the wave are obtained by triadic interactions with the most energetic features of the jet. In our idealised configuration, the wave is deflected in both directions by a given meander. The spatial periodicity of the meander sets the zonal deflections, giving us clearer insight on the BBPOD of the wave knowing only the extended modes of the jet. This method can be seen as complementary to methods based on the bispectrum to detect triadic interactions in the data.

Here we have an a priori on the nature of the interactions, between a low frequency jet and a fast wave field, allowing to study directly their correlations and to propose a methodology well suited for reconstruction. In a future work, we intend to examine the sensitivity of the results to other jet or wave configurations, and even to richer models. Studying the implications of these methods to reconstruct by correlation the internal wave field by the knowledge of the instantaneous slow motion is also envisaged for data assimilation problematics, either by direct estimation, or for learning an observation operator associated with a Galerkin model defined in the BBPOD space.

Figure 1 :

 1 Figure 1: A snapshot taken from the RSW simulation. On the left the vorticity field and on the right the internal wave SSH anomaly h ω , extracted by bandpass filtering centered around the wave frequency ω = 2.

Figure 2 :

 2 Figure 2: Real part of the first three energy-scaled BBPOD modes λ j ψ j associated with u at the tidal frequency. Sponge regions are not shown.

Figure 3 :

 3 Figure 3: Cumulative energy N j=1 λ j (ω)/ +∞ j=1 λ j (ω)[%] contained in the reconstruction of the incoherent field by BBPOD modes.

Figure 4

 4 Figure 4 depicts the direct EPOD modes of the balanced motion, correlated with the BBPOD modes of the wave. The modes are weighted as explained before. The first extended mode is

Figure 4 :

 4 Figure 4: On the left column are presented the real part of the first three BBPOD modes of u at the tidal frequency. On the last column are the corresponding direct EPOD modes of u. Modes are scaled as follows: λ j ψ j .

Figure 5 :

 5 Figure 5: On the first column are presented the first three BBPOD modes of u of the jet. On the last column are the associated inverse EPOD modes of u. Modes are scaled as follows :λ j ψ j .

Table 1 :

 1 Normalized energy of the first 5 modes shown in percentage. Computed in Ω and in Ω N . Each eigenvalue is divided by the total variance in Ω at the tidal frequency j λ Ω j (f ).

	Mode	Ω	Ω N
	1	57.35	29.04
	2	13.81	12.98
	3	9.86	9.51
	4	3.86	3.77
	5	2.34	2.13
	Sum of Incoherence	42.64	37.91