minimum number of sequences for a flank position of 55%. The phylogenetic tree was recon-
structed using the neighbour joining method implemented in the BioNJ program [82] with
N = 100 bootstrapping. The resulting phylogenetic tree was exported as a Newick file and
edited in Evolview [83]. The genomic context around celf1 and celf2 was analysed using the
Genomicus website [84]. A few additional genomes were analysed manually by blasting
(tblastn) some fish NCBI genomes with the spotted gar Celf protein and by manually extract-
ing the NCBI gene annotation around these corresponding celf loci.

Locomotor activity

Locomotor activity was determined at hatching stage (stage 39, 9 dpf) either under dark condi-
tions or with 10% light (approximately equal to 100 lux) by measuring the total distance swum
during a period of 5 minutes (described in [85]) for 12 hatchlings in each condition. Briefly,
the larvae were transferred into 12-well plates containing 1 mL Danieau's solution, and the
plate was placed in a Zebrabox equipped with a video camera, infrared light, and filter (View-
Point Life Sciences, Lyon, France) and the ZebraLab Videotrack software (ViewPoint Life Sci-
ences) for tracking. Following a 5-minute-long habituation period with the same light
conditions as for the experimental period, the locomotor activity was recorded. The detection
threshold was 11, the inactive/small threshold was 0.5 cm/s, and the small/large threshold was
1.0 cm/s. The total distance swum for each individual is equal to the sum of distances reached
during inactivity, small and large movements. *p*-Values were calculated using two-tailed
unpaired *t* test with Excel 14.4.8 (Microsoft), and *p* < 0.05 was considered significant. Bars
and error bars indicate mean ± standard deviation.

EMSA

Binding assays were carried out as previously described [38,40]. In detail, 5'-UGGUUCACGU
CUGCUGCAGGUUCUCUGACUCU-3' for the native D3U-box, 5'-UGGUUCACGUCUGCU
GCAGCUCUCUGACUCU-3' for the Olive fruit fly box (Off-box), 5'-UGGUUCACGUUCU
UCAAGACGCUCUGACUCU-3' for the D3U-scrambled box (S-D3U-box), and 5'-UGGUU
CACGUUCUUCACAGGUUCUCUGACUCU-3' for the D3U-minus CUG box (CUGminus-
D3U-box) were synthesized and end-labelled. For radioactive labelling, 50 pmol of the duplex
5' termini were used, together with 50 pmol of gamma-[³²P]-ATP and 20 units of T4 PNK in
1x-adjusted T4 PNK buffer, and were incubated for 20 minutes at 37 °C. Unincorporated
nucleotides were removed through a Sephadex G-50 spin column. For producing Ol-lrpprc,
Ol-cug-bp1, and Ol-cug-bp2 proteins, pCS2::OL-LRPPRC, pCS2::OL-CUG-BP1 or pCS2::
OL-CUG-BP2 plasmids were linearized using KpnI and then in vitro transcribed using mMes-
sage mMachine kit (Ambion). Finally, Ol-bsf, Ol-cug-bp1, and Ol-cug-bp2 proteins were in
vitro translated using Ambion's Retic Lysate Kit from the previously in vitro-transcribed Ol-
lrpprc, Ol-cug-bp1, and Ol-cug-bp2 capped RNAs. DNA binding reaction contained 10
mMTris-HCl (pH 7.9), 100 mM KCl, 10% glycerol, 5 mM MgCl₂, 1 µg torula rRNA, 0.075%
Triton X-100, 1 mM DTT, 1 µg BSA, 0.5 ng radiolabelled probe, and 2 or 4 µL in vitro transla-
tion mix in a total volume of 20 µL. The amount of 1/10 volume heparin (50 mg/mL) was
added just before loading the binding reaction. For control, reticulocyte lysate alone together
with radiolabelled duplex probe was used and did not result in any shift. Binding reactions
were performed on ice for 10 minutes, and complexes were resolved on a 5% native acrylamide
(37.5:1) gel in 0.5 X TBE and then directly subjected to autoradiography. Of note, due to differ-
ent exposure times—which can be appreciated through the different intensities of the free
probes in the respective figures—autoradiographies should not be compared between each
other.
## Establishment of transgenic fluorescent reporter lines and imaging

For a dynamic and in vivo visualization of endogenous Ol-bsf expression, a transgenic fluorescent reporter line was created. The *Ol-bsf* upstream promoter region of the *Ol-bsf* gene (732 bp up to the next upstream gene) was cloned (BamH1 sites) in front of the mCherry ORF of a meganuclease plasmid (BamH1-BSFp-Fw: AAAGGATCCAGTGTGAGTTCTATCAAGCC TGG; BamH1-BSFp-Rv: AAAGGATCCTTCTGTAGCTGCGTAGAGGAAGATC). For the generation of a stable transgenic line, the meganuclease protocol was used [86]. Briefly, approximately 10 to 15 pg of total vector DNA in a volume of 500 pL injection solution containing *I-SceI* meganuclease was injected into the cytoplasm of 1-cell-staged medaka embryos (Carbio strain). Adult F0 fish were mated to each other, and the offspring were tested for the presence of the transgene by fluorescence check. Siblings from positive F1 fish were raised to adulthood and tested again for fluorescence. Tg[vasa:GFP] and Tg[dmrt1bY:GFP] transgenic lines were described earlier [11,29,48,59].

## Visualization of PGCs

For PGC visualization and counting, the GFP-nos1 3' UTR construct that includes the mmGFP5 ORF cloned upstream of the 3' UTR of the zebrafish *nanos1* gene [5,87] was injected at 1 cell stage. N = 34 and 28 for wild-type and MO-bsf-injected embryos, respectively. For imaging, embryos, hatchlings or tissues were mounted with 1.2% low melting temperature agarose.

## In vivo expression regulation analyses and real-time PCR

For producing Ol-bsf, Ol-cug-bp1, and Ol-cug-bp2 capped mRNAs, pCS2::OL-BSF, pCS2:: OL-CUG-BP1, or pCS2::OL-CUG-BP2 plasmids were linearized using *KpnI* and then in vitro transcribed using mMessage mMachine kit (Ambion). For overexpression, 1 nanolitre was injected into the cytoplasm of 1-cell-stage medaka embryos. For mRNA stability assays, equimolar amounts of control and D3U-box-containing RNAs were injected.

Total RNAs were extracted from fish tissues or embryos using the TRIZOL reagent (Invitrogen) according to the supplier's recommendation. After DNase treatment, reverse transcription was performed with 2 µg total RNA using RevertAid First Strand Synthesis kit (Fermentas) and random primers. Real-time quantitative PCR was carried out with SYBR Green reagents, and amplifications were detected with an i-Cycler (Biorad). All results are averages of at least 2 independent reverse transcription reactions. Error bars represent the standard deviation of the mean. Relative expression levels (according to the equation 2-DeltaCT) were calculated after correction of expression of elongation factor 1 alpha (ef1alpha).

## Morpholino injections

For Ol-bsf knockdown experiments, embryos were injected with a splice morpholino: 5'-TTGATGACTGGCCTGCCAACCTGTC-3' targeting the 3' splice junction of *Ol-bsf* exon 2 (see S3 Fig). The most efficient dose (4 mg/mL) was experimentally determined and the specificity of the targeting confirmed in control experiments (see S3 Fig).

## Microarray

Total RNAs were extracted from adult medaka gonads using the Tri-reagent (Sigma-Aldrich) according to the supplier's recommendation. The total RNA yield was estimated using a NanoDrop ND-1000 spectrophotometer (Labtech, Palaiseau, France), and RNA integrity was checked by means of an Agilent Bioanalyzer (Agilent Technologies, Massy, France). Medaka
gene expression profiling was conducted using an Agilent 8x60K high-density oligonucleotide microarray (GEO platform GPL24100). Labelling and hybridization steps were performed following the 'One-Color Microarray-Based Gene Expression Analysis (Low Input Quick Amp labelling)' Agilent protocol. Briefly, for each sample, 150 ng of total RNA was amplified and labelled using Cy3-CTP. Yield (>825 ng cRNA) and specific activity (>6 pmol of Cy3 per µg of cRNA) of Cy3-cRNA produced were checked with the Nanodrop. The amount of 600 ng of Cy3-cRNA was fragmented and hybridized on a sub-array. Hybridization was carried out for 17 hours at 65 °C in a rotating hybridization oven prior to washing and scanning with an Agilent Scanner (Agilent DNA Microarray Scanner, Agilent Technologies, Massy, France) using the standard parameters for a gene expression 8x60K oligo-array (3 µm and 20 bits). Data were then obtained with the Agilent Feature Extraction software (10.7.3.1) according to the appropriate GE protocol (GE1_107_Sep09). The arrays were normalized (scale normalization) and log-transformed using Genespring Software (version 14.5). A *t* test analysis (*p* < 0.05) was employed to determine the genes that were the most differentially expressed between the 2 conditions. Microarray data sets have been deposited to the GEO-NCBI with the accession number GSE 104726. GO was performed using the panther program (http://geneontology.org/).

## CRISPR/Cas9 genome editing

Identification of CRISPR/Cas9 target sites and design of oligonucleotides of were performed by the use of the ZiFiT software (http://zifit.partners.org/ZiFiT/Disclaimer.aspx). For preparation of sgRNAs, the DR274 plasmid (Addgene number 42250) was first linearized with BsaI, electrophoresed in a 2% agarose gel, and purified. Pairs of complementary oligonucleotides were annealed (40 mM Tris-HCl [pH 8.0], 20 mM MgCl₂, and 50 mM NaCl buffer) by heating at 95 °C for 2 minutes and then cooled down slowly to 25 °C within 1 hour. The double-stranded oligonucleotides were then ligated into the linearized pDR274 vector. Different sgRNAs were designed to target several sites within the Ol-BSF gene (see S6 Fig) in order to create deletions. After linearization with DraI and NotI, respectively, pDR274 and pCS2-nCas9n plasmids were used for generating either sgRNAs or Cas9 RNAs. The synthesized RNAs were then injected into 1-cell-staged embryos at the following concentrations: 25 ng/µL for each sgRNAs and 100 ng/µL for the Cas9 mRNA. CRISPR-positive fish were then screened for mutations using PCR primers flanking the site of deletion (see S6 Fig). The inferred mutant protein is presented in S6C Fig. Mutant fish used in this study have been outcrossed for at least 5 generations.

## Histology

Gonads were fixed for 48 hours in Bouin-Holland fluid and then dehydrated serially in aqueous 70% and 95% ethanol, ethanol/butanol (5:95), and butanol. Tissues were embedded in paraffin, and 5-µm mid-sagittal gonad sections were stained with Regaud's haematoxylin and haematoxylin-eosin-safran.

## Supporting information

**S1 Fig. Presence of the D3U-box in different UTRs and coding sequence of gonadal genes.**
Presence of the D3U-box (matrix-scan) was evaluated within the transcriptome of the medaka fish. D3U-box, *dmrt1* 3' UTR box.
(DOCX)
**S2 Fig. Phylogeny and synteny analysis of Ol-bsf and Ol-cug-bp ohnologs.** (A) Circular cladogram representation of the phylogenetic tree of Lrpprc proteins in jawed vertebrates (gnatosthomes). This phylogeny shows that *lrpprc* genes were retained as single copies in most jawed vertebrates even following whole genome duplications (red stars), e.g., the teleost-specific duplication or the salmonid-specific duplication. *Lrpprc* is, however, present in duplicated copies in the Cyprininae (tree branches in red). Bootstraps (*N* = 100) values are indicated in each tree node when judged significant (>0.75). Tree branches are depicted in blue for lobefinal vertebrates and cartilaginous fish and in black for teleosts with the exception of Cyprininae in red. (B) Gene evolution of *celf2* genes in some teleosts. The phylogeny on the left is a dendogram representation of *celf2* gene phylogeny in teleosts given as an indication as only a few nodes are supported by good bootstraps' values (*N* = 100, mentioned in each tree nodes when judged significant, i.e., >0.7). The teleost fish whole genomic duplication (3R) is indicated by a red star. The left part of the figure is a representation of the evolution of the genomic context around the *celf2* gene. After the 3R whole genome duplication, *celf2*—which is a single copy gene on the Chr 8 of the spotted gar genome—was duplicated in two 3R ohnologs, *celf2a* and *celf2b*, that were not retained as 2 copies in all teleost fish. The genomic context of the *celf2a* and *celf2b* paralogous regions clearly indicates a partition of the ancestral region found in spotted gar. The *celf2a* gene was retained in all species investigated, but the *celf2b* gene seems to have been lost in Otophysi or at least in *D. rerio* (Cypriniformes), *Astyanax mexicanus* (Characiformes), and *Ictalurus punctatus* (Siluriformes). *celf2*, CUGBP Elav-like family member 2; Chr 8, Chromosome 8; Lrpprc, leucine rich pentatricopeptide repeat containing; Ol-BSF, *Oryzias latipes* Bicoid Stability Factor; Ol-CUG-BP, *Oryzias latipes* CUG-binding protein. (TIF)

**S3 Fig. Analysis of morpholino efficiency and level of Ol-bsf down-regulation.** For in vivo transient down-regulation of Ol-bsf, a splice morpholino was designed to encompass the splice junction between exon 2 and intron 2 of the *Ol-bsf* gene in order to induce aberrant splicing and frame shit of the ORF. To show to what extend the splicing/activity of *Ol-bsf* was impacted, RT-PCR using exons 1, 2, and 3 spanning primers together with cDNAs from different stages of morpholino-injected embryos was achieved. E2, exon 2; i2, intron 2; Ol-BSF, *Oryzias latipes* Bicoid Stability Factor; RT-PCR, Reverse Transcription-Polymerase Chain Reaction.

(TIF)

**S4 Fig. Real-time PCR quantification of Ol-cug-bp1, Ol-cug-bp2, and Ol-bsf during embryogenesis and in adult tissues.** (A and C) During embryonic development, both Ol-cug-bp ohnologs are expressed in a complementary manner. Being likely maternally deposited the expression of Ol-cug-bp1 rapidly decreases after mid-blastula transition (stage 10) to remain virtually off up to hatching stage. On the other hand, low expression of Ol-cug-bp2 is detected until stage 25 and rapidly increases by stage 33. (B and D) In adult tissues, both Ol-cug-bp ohnologs are expressed in brain, muscles, and gonads; ol-cug-bp2 is additionally expressed in eyes and skin. Both ohnologs are higher expressed in male gonads than in female gonads. (E and F) In adult tissues, Ol-bsf is ubiquitously present although higher expression is observed in gonads of both sexes. Underlying data for (A to F) can be found in S1 Data.

(TIF)

**S5 Fig. Lrrprc and celf1, but not celf2, are expressed in mouse embryonic gonads.** (A to H) ISHs on sagittal sections of 14.5 dpc mouse embryos showed expression of Lrrprc (A to D) and celf1 (E to H) most likely in germ cells within testis cords (B and F) and germ cells within the ovary (D and H). In contrast, no celf2 expression was detected in developing gonads (I–L).
However, celf2 expression was detected in other tissues, such as part of the brain and dorsal root ganglia. Scale bars: 1 mm for A, C, E, G, I, and K; 10 mm for B, D, F, H, J, and L. (M–R) Immunofluorescent detection of LRPPRC (M, N, P, Q) and DDX4/VASA (O, R) in adult mouse testes (M–O) and ovaries (P–R). In adult testes, lrpprc is expressed in one subpopulation of germ cells; compared lrpprc staining on (M) and (N) with vasa staining on (O) where most of the germ cells (except some spermatogonia) remain stained by vasa. According to the position of lrpprc-positive cells (arrowheads in M and N) in the seminiferous tubule (not basal and below round spermatids) and to the fact that lrpprc-positive germ cells are those with the largest nucleus, lrpprc-positive cells seem to be spermatocytes at the pachytene stage. In adult ovaries (P–R), lrpprc is mainly expressed into the oocytes of primordial, primary and young secondary follicles (see arrows on [P] and [Q]). LRPPRC staining disappears from the oocyte of secondary follicles that are clearly stained for vasa in (R) (compared stars in [Q] and [R]). Scale bars: 200 µm for M to R. dph, days post hatching; ISH, in situ hybridisation; LRPPRC, leucine rich pentatricopeptide repeat containing.

(TIF)

**S6 Fig. Generation *Ol-bsf* knockout medaka fish after genome editing by CRISPR/Cas9 method.** (A) Several guide RNA were designed in order to target different locations on the *Ol bsf* gene (targets 1, 2, and 6). (B) After injection of different combinations of guide RNAs together with the Cas9 mRNA, putative edited fish were subjected to RT-PCR using primer sets flanking the cutting sites. Lines displaying deletions within the *Ol bsf* gene (red stars) were kept for further investigations. (C) Deletions obtained within the *Ol bsf* gene (left panel) result in a truncated translated Ol-bsf protein (right panel). CRISPR-Cas9, clustered regularly inter-spaced short palindromic repeats/CRISPR-associated protein 9; Ol-BSF, *Oryzias latipes* Bicoid Stability Factor; RT-PCR, Reverse Transcription- Polymerase Chain Reaction.

(TIF)

**S7 Fig. Locomotor activity of *Ol-bsf* mutant fish versus wild type.** Locomotor activity (*Ol bsf* mutants versus wild type) was determined at hatching stage (stage 39, 9 dpf) either under dark conditions (A) or with 10% light (approximately equal to 100 lux [panel B]) by measuring the total distance swum during a period of 5 minutes. (C) The total distance swum for each individual is equal to the sum of distances reached during inactivity, small and large movements. Bars and error bars indicate mean ± standard deviation. N = 12 for each condition. Underlying data for (C) can be found in S1 Data. dpf, days post fertilization; Ol-BSF, *Oryzias latipes* Bicoid Stability Factor.

(TIF)

**S8 Fig. Ovarian phenotypes of the *Ol-bsf* mutant fish.** Morphological inspection of heterozygote mutant ovaries discloses a significant accumulation of small sized-oocytes compared to wild type (A) and (A1 to A5 for details and statistical analyses). (A) Overall size distribution of the oocytes in 9 wild-type and 9 Ol-BSF(+/-) adult ovaries. Each gonad (testes or ovaries) was sectioned through the mid-sagittal plan (see also Materials and methods). Underlying data for (A) can be found in S1 Data. Ol-BSF, *Oryzias latipes* Bicoid Stability Factor.

(TIF)

**S9 Fig. Ovarian phenotypes of the *Ol-BSF* mutant fish.** (A to R) Mid-sagittal sections of the ovaries utilized for counting the oocytes in S8 Fig. Each gonad (testes or ovaries) was sectioned through the mid-sagittal plan (see also Materials and methods). Underlying data for (A) can be found in S1 Data. Ol-BSF, *Oryzias latipes* Bicoid Stability Factor.

(TIF)
**S10 Fig. Testicular phenotypes of the Ol-bsf mutant fish.** (A–J compared to K–T) Heterozygote mutant testes (A–J) exhibit a decreased number of spermatogonia with accumulation of type 2 spermatocytes, spermatids, and sperm within the most external layers of the seminiferous epithelium (arrowheads) compared to wild-type testes (K–T). Either 10 different wild-type (A–J) or *Ol-bsf*-deficient (K–T) testes were analysed. Mid-sagittal gonad sections were stained with haematoxylin–eosin–safran. Each gonad (testes or ovaries) was sectioned through the mid-sagittal plan (see also *Material and methods*). Ol-BSF, *Oryzias latipes* Bicoid Stability Factor.
(TIF)

**S11 Fig. Fertility test.** Egg numbers and fertilization rates were recorded over a period of 9 days for the following crosses: (A) male *Ol-bsf*(-/+) × female *Ol-bsf*(-/+); (B) male *Ol-bsf*(-/+) × female wild type; (C) male wild type × female *Ol-bsf*(-/+). Underlying data for (A to C) can be found in S1 Data. BSF, bicoid stability factor.
(TIF)

**S12 Fig. Microarray data and mitochondrial gene quantification.** (A) Microarray. Adult testes of either *bsf*⁺/⁻ or wild-type animals were subjected to microarray (see *Materials and methods*). GO term analysis reveals that in mutant testes partially depleted for the *bsf* gene, rRNA processing is particularly affected. Ol-bsf and Ol-cug-bp2 are down- and up-regulated, respectively, in mutant animals compared to wild type. Of note, and in accordance with the literature, a significant proportion (10.1%) of the down-regulated genes code for proteins localized in the mitochondria. Finally, supporting our observations that lowering *ol-bsf* transcription (morpholino injection in Fig 7E) resulted in up-modulation of germ cell number and that mutant gonads presented an increase of germ cells committing to gametogenesis (Fig 6), our microarray analysis reveals a general up-regulation of genes involved in germ cell proliferation or differentiation. (B) RNA levels of different mitochondrial genes (Cox1, Cox2, ND1, ND5, and CytB) were quantified by real-time PCR after BSF-morpholino injections and compared to wild type. Most of the mitochondrial genes are down-regulated when the level of Ol-bsf decreases. (C) Modulation of RNA levels of the cyp19a1 (aromatase) gene after overexpression of Ol-cugbp1 or Ol-cug-bp2. Underlying data for (B and C) can be found in S1 Data. GO, gene ontology.
(TIFF)

**S1 Table. Gene evolution of *cugbp elav-like family member 2* genes in some teleosts. *celf2, cugbp elav-like family member 2*.
(DOCX)**

**S1 Data. Underlying data.**
(XLSX)

**Acknowledgments**

Acknowledgements

The authors thank the fish facilities at the University of Würzburg and at the INRA Rennes (C. Duret, F. Borel, and A. Patinote).

**Author Contributions**

**Conceptualization:** Amaury Herpin, Manfred Schartl.

**Formal analysis:** Amaury Herpin, Clara Gobé, Aurélie Le Cam, Yann Guiguen.

**Funding acquisition:** Amaury Herpin.
**Investigation:** Amaury Herpin, Cornelia Schmidt, Clara Gobé, Martina Regensburger, Aurélie Le Cam, Jérome Montfort, Mateus C. Adolfi, Christina Lillesaar, Dagmar Wilhelm, Michael Kraeusling, Brigitte Mourot, Béatrice Porcon, Laurence Ettwiller, Dirk Dolle, Yann Guiguen, Manfred Schartl.

**Methodology:** Amaury Herpin, Christina Lillesaar, Maëlle Pannetier, Eric Pailhoux, Laurence Ettwiller, Dirk Dolle.

**Project administration:** Amaury Herpin.

**Software:** Susanne Kneitz, Laurence Ettwiller, Dirk Dolle.

**Supervision:** Amaury Herpin, Manfred Schartl.

**Validation:** Amaury Herpin.

**Writing – original draft:** Amaury Herpin, Manfred Schartl.

**Writing – review & editing:** Amaury Herpin, Manfred Schartl.

## References

1. Brennan J, Capel B. One tissue, two fates: molecular genetic events that underlie testis versus ovary development. Nat Rev Genet. 2004; 5: 509–521. https://doi.org/10.1038/nrg1381 PMID: 15211353

2. Devlin RH, Nagahama Y. Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture. 2002; 208: 191–364. https://doi.org/10.1016/S0044-8486(02)00057-1

3. Lin Y-T, Capel B. Cell fate commitment during mammalian sex determination. Curr Opin Genet Dev. 2015; 32: 144–152. https://doi.org/10.1016/j.gde.2015.03.003 PMID: 25841206

4. Rios-Rojas C, Bowles J, Koopman P. On the role of germ cells in mammalian gonad development: quiet passengers or back-seat drivers? Reprod Camb Engl. 2015; 149: R181–191. https://doi.org/10.1530/REP-14-0663 PMID: 25628441

5. Herpin A, Fischer P, Liedtke D, Kluever N, Neuner C, Raz E, et al. Sequential SDF1α and b-induced mobility guides Medaka PGC migration. Dev Biol. 2008; 320: 319–327. https://doi.org/10.1016/j.ydbio.2008.03.030 PMID: 18440502

6. Kurokawa H, Saito D, Nakamura S, Katoh-Fukui Y, Ohta K, Baba T, et al. Germ cells are essential for sexual dimorphism in the medaka gonad. Proc Natl Acad Sci U S A. 2007; 104: 16958–16963. https://doi.org/10.1073/pnas.0609932104 PMID: 17940041

7. Marlow F. Primordial Germ Cell Specification and Migration. F1000Research. 2015; 4. https://doi.org/10.12688/f1000research.6995.1 PMID: 26918157

8. Nakamura S, Kobayashi D, Aoki Y, Yokoi H, Ebe Y, Wittbrodt J, et al. Identification and lineage tracing of two populations of somatic gonadal precursors in medaka embryos. Dev Biol. 2006; 295: 678–688. https://doi.org/10.1016/j.ydbio.2006.03.052 PMID: 16682019

9. Paksa A, Raz E. Zebrafish germ cells: motility and guided migration. Curr Opin Cell Biol. 2015; 36: 80–85. https://doi.org/10.1016/j.ceb.2015.07.007 PMID: 26232877

10. Johnson AD, Richardson E, Bachvarova RF, Crother BI. Evolution of the germ line-soma relationship in vertebrate embryos. Reprod Camb Engl. 2011; 141: 291–300. https://doi.org/10.1530/REP-10-0474 PMID: 21228047

11. Tanaka M, Kinoshita M, Kobayashi D, Nagahama Y. Establishment of medaka (Oryzias latipes) transgenic lines with the expression of green fluorescent protein fluorescence exclusively in germ cells: a useful model to monitor germ cells in a live vertebrate. Proc Natl Acad Sci U S A. 2001; 98: 2544–2549. https://doi.org/10.1073/pnas.041315498 PMID: 11226275

12. Herpin A, Schartl M. Dtmt1 genes at the crossroads: a widespread and central class of sexual development factors in fish. FEBS J. 2011; 278: 1010–1019. https://doi.org/10.1111/j.1742-4658.2011.08030.x PMID: 21281449

13. Herpin A, Schartl M. Dtmt1 genes at the crossroads: a widespread and central class of sexual development factors in fish. FEBS J. 2011; 278: 1010–1019. https://doi.org/10.1111/j.1742-4658.2011.08030.x PMID: 21281449

14. Matson CK, Zarkower D. Sex and the singular DM domain: insights into sexual regulation, evolution and plasticity. Nat Rev Genet. 2012; 13: 163–174. https://doi.org/10.1038/nrg3161 PMID: 22310892
15. Raymond CS, Murphy MW, O'Sullivan MG, Bardwell VJ, Zarkower D. Dmrt1, a gene related to worm and fly sexual regulators, is required for mammalian testis differentiation. *Genes Dev.* 2000; 14: 2587–2595. PMID: 11040213

16. Chong T, Collins JJ, Brubacher JL, Zarkower D, Newmark PA. A sex-specific transcription factor controls male identity in a simultaneous hermaphrodite. *Nat Commun.* 2013; 4: 1814. https://doi.org/10.1038/ncomms2811 PMID: 23652002

17. Kato Y, Kobayashi K, Watanabe H, Iguchi T. Environmental sex determination in the branchiopod crustacean Daphnia magna: deep conservation of a Doublesex gene in the sex-determining pathway. *PLoS Genet.* 2011; 7: e1001345. https://doi.org/10.1371/journal.pgen.1001345 PMID: 21455482

18. Yoshimoto S, Okada E, Umemoto H, Tamura K, Uno Y, Nishida-Umehara C, et al. A W-linked DM-domain gene, DM-W, participates in primary ovary development in Xenopus laevis. *Proc Natl Acad Sci U S A.* 2008; 105: 2469–2474. https://doi.org/10.1073/pnas.0712244105 PMID: 18268317

19. Chen S, Zhang G, Shao C, Huang Q, Liu G, Zhang P, et al. Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. *Nat Genet.* 2014; 46: 253–260. https://doi.org/10.1038/ng.2890 PMID: 24487278

20. Smith CA, Roeszler KN, Ohnesorg T, Cummins DM, Farlie PG, Doran TJ, et al. The avian Z-linked gene DMRT1 is required for male sex determination in the chicken. *Nature*. 2009; 461: 267–271. https://doi.org/10.1038/nature08298 PMID: 19710650

21. Matsuda M, Nagahama Y, Shinomiya A, Sato T, Matsuda C, Kobayashi T, et al. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. *Nature*. 2002; 417: 559–563. https://doi.org/10.1038/nature751 PMID: 12037570

22. Nanda I, Kondo M, Hornung U, Asakawa S, Winkler C, Shimizu A, et al. A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes. *Proc Natl Acad Sci U S A.* 2002; 99: 11778–11783. https://doi.org/10.1073/pnas.182314699 PMID: 12193652

23. Ledig S, Hiort O, Wünsch L, Wieacker P. Partial deletion of DMRT1 causes 46,XY ovotesticular disorder of sexual development. *Eur J Endocrinol.* 2012; 167: 119–124. https://doi.org/10.1530/EJE-12-0136 PMID: 22573722

24. Raymond CS, Parker ED, Kettlewell JR, Brown LG, Page DC, Kusz K, et al. A region of human chromosome 9p required for testis development contains two genes related to known sexual regulators. *Hum Mol Genet.* 1999; 8: 989–996. PMID: 10332030

25. Kim S, Bardwell VJ, Zarkower D. Cell type-autonomous and non-autonomous requirements for Dmrt1 in postnatal testis differentiation. *Dev Biol.* 2007; 307: 314–327. https://doi.org/10.1016/j.ydbio.2007.04.046 PMID: 17540358

26. Matson CK, Murphy MW, Griswold MD, Yoshida S, Bardwell VJ, Zarkower D. The mammalian double-sex homolog DMRT1 is a transcriptional gatekeeper that controls the mitosis versus meiosis decision in male germ cells. *Dev Cell.* 2010; 19: 612–624. https://doi.org/10.1016/j.devcel.2010.09.010 PMID: 20951351

27. Matson CK, Murphy MW, Sarver AL, Griswold MD, Bardwell VJ, Zarkower D. DMRT1 prevents female reprogramming in the postnatal mammalian testis. *Nature*. 2011; 476: 101–104. https://doi.org/10.1038/nature10239 PMID: 21775990

28. Herpin A, Schartl M. Plasticity of gene-regulatory networks controlling sex determination: of masters, slaves, usual suspects, newcomers, and usurpators. *EMBO Rep.* 2015; 16: 1260–1274. https://doi.org/10.15252/embr.201540667 PMID: 26358957

29. Herpin A, Schindler D, Kraiss A, Hornung U, Winkler C, Schartl M. Inhibition of primordial germ cell proliferation by the medaka male determining gene Dmrt I bY. *BMC Dev Biol.* 2007; 7: 99. https://doi.org/10.1186/1471-213X-7-99 PMID: 17760954

30. Krentz AD, Murphy MW, Kim S, Cook MS, Capel B, Zhu R, et al. The DM domain protein DMRT1 is a dose-sensitive regulator of fetal germ cell proliferation and pluripotency. *Proc Natl Acad Sci U S A.* 2009; 106: 22323–22328. https://doi.org/10.1073/pnas.0905431106 PMID: 20007774

31. Shao C, Li Q, Chen S, Zhang P, Lian J, Hu Q, et al. Epigenetic modification and inheritance in sexual reversal of fish. *Genome Res.* 2014; 24: 604–615. https://doi.org/10.1101/gr.162172.113 PMID: 24487721

32. Wen AY, You F, Sun P, Li J, Xu DD, Wu ZH, et al. CpG methylation of dmrt1 and cyp19a promoters in relation to their sexual dimorphic expression in the Japanese flounder *Paralichthys olivaceus*. *J Fish Biol.* 2014; 84: 193–205. https://doi.org/10.1111/jfb.12277 PMID: 24372528

33. Ge C, Ye J, Zhang H, Zhang Y, Sun W, Sang Y, et al. Dmrt1 induces the male pathway in a turtle species with temperature-dependent sex determination. *Dev Camb Engl.* 2017; 144: 2222–2233. https://doi.org/10.1242/dev.152033 PMID: 28506988
34. Lei N, Heckert LL. Sp1 and Egr1 regulate transcription of the Dmrt1 gene in Sertoli cells. Biol Reprod. 2002; 66: 675-684. PMID: 11870074

35. Lei N, Heckert LL. Gata4 regulates testis expression of Dmrt1. Mol Cell Biol. 2004; 24: 377-388. https://doi.org/10.1128/MCB.24.1.377-388.2004 PMID: 14673170

36. Cui J, Shen X, Zhao H, Nagahama Y. Genome-wide analysis of Sox genes in Medaka (Oryzias latipes) and their expression pattern in embryonic development. Cytogenet Genome Res. 2011; 134: 283-294. https://doi.org/10.1159/000329480 PMID: 21720154

37. Matsuda M, Shinomiya A, Kinoshita M, Suzuki A, Kobayashi T, Paul-Prasanth B, et al. DMY gene induces male development in genetically female (XX) medaka fish. Proc Natl Acad Sci U S A. 2007; 104: 3865-3870. https://doi.org/10.1073/pnas.0611707104 PMID: 17360444

38. Herpin A, Braasch I, Kraeussling M, Schmidt C, Thoma EC, Nakamura S, et al. Transcriptional rewiring of the sex determining dmrt1 gene duplicate by transposable elements. PLoS Genet. 2010; 6: e1000844. https://doi.org/10.1371/journal.pgen.1000844 PMID: 20169179

39. Schartl M, Schories S, Wakamatsu Y, Nagao Y, Hashimoto H, Bertin C, et al. Sox5 is involved in germ-cell regulation and sex determination in medaka following co-option of nested transposable elements. BMC Biol. 2018; 16: 16. https://doi.org/10.1186/s12915-018-0485-8 PMID: 29378592

40. Herpin A, Nakamura S, Wagner TU, Tanaka M, Schartl M. A highly conserved cis-regulatory motif directs differential gonadal synexpression of Dmrt1 transcripts during gonad development. Nucleic Acids Res. 2009; 37: 1510-1520. https://doi.org/10.1093/nar/gkn1065 PMID: 19139075

41. Mancebo R, Zhou X, Shillinglaw W, Henzel W, Macdonald PM. BSF binds specifically to the bicoid mRNA 3' untranslated region and contributes to stabilization of bicoid mRNA. Mol Cell Biol. 2001; 21: 3462-3471. https://doi.org/10.1128/MCB.21.10.3462-3471.2001 PMID: 11313472

42. De Renzis S, Elemento O, Tavazoie S, Wieschaus EF. Unmasking activation of the zygotic genome using chromosomal deletions in the Drosophila embryo. PLoS Biol. 2007; 5: e117. https://doi.org/10.1371/journal.pbio.0050117 PMID: 17456005

43. Paillard L, Omilli F, Legagneux V, Bassez T, Maniey D, Osborne HB. EDEN and EDEN-BP, a cis element and an associated factor that mediate sequence-specific mRNA deadenylation in Xenopus embryos. EMBO J. 1998; 17: 278-287. https://doi.org/10.1093/emboj/17.1.278 PMID: 9427761

44. Philips AV, Timchenko LT, Cooper TA. Disruption of splicing regulated by a CUG-binding protein in myotonic dystrophy. Science. 1998; 280: 737-741. PMID: 9563950

45. Timchenko LT, Miller JW, Timchenko NA, DeVore DR, Datar KV, Lin L, et al. Identification of a (CUG)n triplet repeat RNA-binding protein and its expression in myotonic dystrophy. Nucleic Acids Res. 1996; 24: 4407-4414. PMID: 8948631

46. Paillard L, Legagneux V, Maniey D, Osborne HB. c-Jun ARE targets mRNA deadenylation by an EDEN-BP (embryo deadenylation element-binding protein)-dependent pathway. J Biol Chem. 2002; 277: 3232-3235. https://doi.org/10.1074/jbc.M109362200 PMID: 11707455

47. Moraes KCM, Wilusz CJ, Wilusz J. CUG-BP binds to RNA substrates and recruits PARN deadenylase. RNA N Y N. 2006; 12: 1084-1091. https://doi.org/10.1261/rna.59606 PMID: 16601207

48. Hornung U, Herpin A, Schartl M. Expression of the male determining gene dmrt1bY and its autosomal coorthologue dmrt1a in medaka. Sex Dev Genet Mol Biol Evol Endocrinol Embryol Pathol Sex Determ Differ. 2007; 1: 197-206. https://doi.org/10.1159/000102108 PMID: 18391530

49. Webster KA, Schach U, Ordaz A, Steinfeld JS, Draper BW, Siegfried KR. Dmrt1 is necessary for male sexual development in zebrafish. Dev Biol. 2017; 422: 33-46. https://doi.org/10.1016/j.ydbio.2016.12.008 PMID: 27940159

50. Shiraishi E, Imazato H, Yamamoto T, Yokoi H, Abe S, Kitano T. Identification of two teleost homologs of the Drosophila sex determination factor, transformer-2 in medaka (Oryzias latipes). Mech Dev. 2004; 121: 991-996. https://doi.org/10.1016/j.mod.2004.04.013 PMID: 15210204

51. Amrein H, Gorman M, Nöthiger R. The sex-determining gene tra-2 of Drosophila encodes a putative RNA binding protein. Cell. 1988; 55: 1025-1035. PMID: 3144434

52. Shiraishi E, Yoshinaga N, Miura T, Yokoi H, Wakamatsu Y, Abe S-I, et al. Mullerian inhibiting substance is required for germ cell proliferation during early gonadal differentiation in medaka (Oryzias latipes). Endocrinology. 2008; 149: 1813-1819. https://doi.org/10.1210/en.2007-1535 PMID: 18162517

53. Weidinger G, Stebler J, Slanchev K, Dumstrei K, Wise C, Lovell-Badge R, et al. dead end, a novel vertebrate germ plasm component, is required for zebrafish primordial germ cell migration and survival. Curr Biol CB. 2003; 13: 1429-1434. PMID: 12932328

54. Gross-Thebing T, Yigit S, Pfeiffer J, Reichman-Fried M, Bandemer J, Ruckert C, et al. The Vertebrate Protein Dead End Maintains Primordial Germ Cell Fate by Inhibiting Somatic Differentiation. Dev Cell. 2017; 43: 704-715.e5. https://doi.org/10.1016/j.devcel.2017.11.019 PMID: 29257950
55. Hartung O, Forbes MM, Marlow FL. Zebrafish vasa is required for germ-cell differentiation and maintenance. Mol Reprod Dev. 2014; 81: 946-961. https://doi.org/10.1002/mrd.22414 PMID: 25257909

56. Siegfried KR, Nüsslein-Volhard C. Germ line control of female sex determination in zebrafish. Dev Biol. 2008; 324: 277-287. https://doi.org/10.1016/j.ydbio.2008.09.025 PMID: 18930041

57. Siegfried KR. In search of determinants: gene expression during gonadal sex differentiation. J Fish Biol. 2010; 76: 1879-1902. https://doi.org/10.1111/j.1095-8649.2010.02594.x PMID: 20557645

58. Nakamura S, Watakabe I, Nishimura T, Toyoda A, Taniguchi Y, Tanaka M. Analysis of medaka sox9 orthologue reveals a conserved role in germ cell maintenance. PLoS ONE. 2012; 7: e29982. https://doi.org/10.1371/journal.pone.0029982 PMID: 22253846

59. Nishimura T, Herpin A, Kimura T, Hara I, Kawasaki T, Nakamura S, et al. Analysis of a novel gene, Sdgc, reveals sex chromosome-dependent differences of medaka germ cells prior to gonad formation. Dev Camb Engl. 2014; 141: 3363-3369. https://doi.org/10.1242/dev.106864 PMID: 25078651

60. Herpin A, Rohr S, Riedel D, Kluever N, Raz E, Schartl M. Specification of primordial germ cells in medaka (Oryzias latipes). BMC Dev Biol. 2007; 7: 3. https://doi.org/10.1186/1471-213X-7-3 PMID: 17217535

61. Herpin A, Rohr S, Riedel D, Kluever N, Raz E, Schartl M. Specification of primordial germ cells in medaka (Oryzias latipes). BMC Dev Biol. 2007; 7: 3. https://doi.org/10.1186/1471-213X-7-3 PMID: 17217535

62. Herpin A, Adolfi MC, Nicol B, Hinzmann M, Schmidt C, Klughammer J, et al. Divergent expression regulation of gonad development genes in medaka shows incomplete conservation of the downstream regulatory network of vertebrate sex determination. Mol Biol Evol. 2013; 30: 2328-2346. https://doi.org/10.1093/molbev/mst130 PMID: 23883523

63. Nakamura S, Kobayashi K, Nishimura T, Higashijima S, Tanaka M. Identification of germline stem cells in the ovary of the teleost medaka. Science. 2010; 328: 1561-1563. https://doi.org/10.1126/science.1185473 PMID: 20488987

64. Uribe MC, Grier HJ, Mejía-Roa V. Comparative testicular structure and spermatogenesis in bony fishes. Spermatogenesis. 2015; 4. https://doi.org/10.4161/21565562.2014.983400 PMID: 26413405

65. Keene JD. RNA regulons: coordination of post-transcriptional events. Nat Rev Genet. 2007; 8: 533-543. https://doi.org/10.1038/nrg2111 PMID: 17572691

66. Moore MJ. From birth to death: the complex lives of eukaryotic mRNAs. Science. 2005; 309: 1514-1518. https://doi.org/10.1126/science.1111443 PMID: 16141059

67. Arvola RM, Weidmann CA, Tanaka Hall TM, Goldstrohm AC. Combinatorial control of messenger RNAs by Pumilio, Nanos and Brain Tumor Proteins. RNA Biol. 2017; 1-12. https://doi.org/10.1080/15476286.2017.1306168 PMID: 28318367

68. Lasko P. Posttranscriptional regulation in Drosophila oocytes and early embryos. Wiley Interdiscip Rev RNA. 2011; 2: 408-416. https://doi.org/10.1002/wrna.70 PMID: 21957026

69. Lehmann R. Germ Plasm Biogenesis—An Oskar-Centric Perspective. Curr Top Dev Biol. 2016; 116: 679-707. https://doi.org/10.1016/bs.ctdb.2015.11.024 PMID: 26970648

70. Herpin A, Adolfi MC, Nicol B, Hinzmann M, Schmidt C, Klughammer J, et al. Divergent expression regulation of gonad development genes in medaka shows incomplete conservation of the downstream regulatory network of vertebrate sex determination. Mol Biol Evol. 2013; 30: 2328-2346. https://doi.org/10.1093/molbev/mst130 PMID: 23883523

71. Sterky FH, Ruzzenente B, Gustafsson CM, Samuelsson T, Larsson N-G. LRPPRC is a mitochondrial matrix protein that is conserved in metazoans. Biochem Biophys Res Commun. 2010; 398: 759-764. https://doi.org/10.1016/j.bbrc.2010.07.019 PMID: 20633537

72. Kress C, Gautier-Courteille C, Osborne HB, Babinet C, Paillard L. Inactivation of CUG-BP1/CELF1 causes growth, viability, and spermatogenesis defects in mice. Mol Cell Biol. 2007; 27: 1146-1157. https://doi.org/10.1128/MCB.01009-06 PMID: 17130239

73. Gohil VM, Nilsson R, Belcher-Timme CA, Luo B, Root DE, Mootha VK. Mitochondrial and nuclear genomic responses to loss of LRPPRC expression. J Biol Chem. 2010; 285: 13742-13747. https://doi.org/10.1074/jbc.M109.098400 PMID: 20220140

74. Sasarman F, Brunel-Guitton C, Antonicka H, Wai T, Shoubridge EA, LSFC Consortium. LRPPRC and SLIRP interact in a ribonucleoprotein complex that regulates posttranscriptional gene expression in mitochondria. Mol Biol Cell. 2010; 21: 1315-1323. https://doi.org/10.1091/mbc.E10-01-0047 PMID: 20200222

75. Bratic A, Wredenberg A, Grönke S, Stewart JB, Mourier A, Ruzzenente B, et al. The bicoid stability factor controls polyadenylation and expression of specific mitochondrial mRNAs in Drosophila melanogaster. PLoS Genet. 2011; 7: e1002324. https://doi.org/10.1371/journal.pgen.1002324 PMID: 22022283
76. Sasarman F, Nishimura T, Antonicka H, Weraarpachai W, Shoubridge EA, LSFC Consortium. Tissue-specific responses to the LRPPRC founder mutation in French Canadian Leigh Syndrome. Hum Mol Genet. 2015; 24: 480-491. https://doi.org/10.1093/hmg/ddu468 PMID: 25214534

77. Boulanger G, Cibois M, Viet J, Fostier A, Deschamps S, Pastezeur S, et al. Hypogonadism Associated with Cyp19a1 (Aromatase) Posttranscriptional Upregulation in Celf1 Knockout Mice. Mol Cell Biol. 2015; 35: 3244-3253. https://doi.org/10.1128/MCB.00074-15 PMID: 26169831

78. Iwamatsu T. Stages of normal development in the medaka Oryzias latipes. Mech Dev. 2004; 121: 605-618. https://doi.org/10.1016/j.mod.2004.03.012 PMID: 15210170

79. Hong Y, Winkler C, Schartl M. Production of medakafish chimeras from a stable embryonic stem cell line. Proc Natl Acad Sci U S A. 1998; 95: 3679-3684. PMID: 9520425

80. Dereeper A, Guignon V, Blanc G, Audic S, Buffet S, Chevenet F, et al. Phylogeny.fr: robust phylogenetic analysis for the non-specialist. Nucleic Acids Res. 2008; 36: W465-W469. https://doi.org/10.1093/nar/gkn180 PMID: 18424797

81. Pasquier J, Cabau C, Nguyen T, Jouanno E, Severac D, Braasch I, et al. Gene evolution and gene expression after whole genome duplication in fish: the PhyloFish database. BMC Genomics. 2016; 17: 368. https://doi.org/10.1186/s12864-016-2709-z PMID: 27189481

82. Gascuel O. BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Mol Biol Evol. 1997; 14: 685-695. https://doi.org/10.1093/oxfordjournals.molbev.a025808 PMID: 9254330

83. Zhang H, Gao S, Lercher MJ, Hu S, Chen W-H. EvolView, an online tool for visualizing, annotating and managing phylogenetic trees. Nucleic Acids Res. 2012; 40: W569-W572. https://doi.org/10.1093/nar/gks576 PMID: 22695796

84. Louis A, Muffato M, Roest Crollius H. Genomicus: five genome browsers for comparative genomics in eukaryota. Nucleic Acids Res. 2013; 41: D700-D705. https://doi.org/10.1093/nar/gks1156 PMID: 23193262

85. Lange M, Norton W, Coolen M, Chaminade M, Merker S, Proft F, et al. The ADHD-linked gene Lphn3.1 controls locomotor activity and impulsivity in zebrafish. Mol Psychiatry. 2012; 17: 855. https://doi.org/10.1038/mp.2012.119 PMID: 22918194

86. Grabher C, Wittbrodt J. Meganuclease and transposon mediated transgenesis in medaka. Genome Biol. 2007; 8 Suppl 1: S10. https://doi.org/10.1186/gb-2007-8-s1-s10 PMID: 18047687

87. Köprunner M, Thisse C, Thisse B, Raz E. A zebrafish nanos-related gene is essential for the development of primordial germ cells. Genes Dev. 2001; 15: 2877-2885. https://doi.org/10.1101/gad.212401 PMID: 11691838
# Reconstruction of the birth of a male sex chromosome present in Atlantic herring

Nima Rafati<sup>a,b,1</sup>, Junfeng Chen<sup>a,1</sup>, Amaury Herpin<sup>c,d,1</sup>, Mats E. Pettersson<sup>a,1</sup>, Fan Han<sup>a</sup>, Chungang Feng<sup>a</sup>, Ola Wallerman<sup>a</sup>, Carl-Johan Rubin<sup>a</sup>, Sandrine Péron<sup>c</sup>, Arianna Cocco<sup>a</sup>, Mårten Larsson<sup>a</sup>, Christian Trötschel<sup>e</sup>, Ansgar Poetsch<sup>e,f,g</sup>, Kai Korsching<sup>h</sup>, Wolfgang Bönigk<sup>h</sup>, Heinz G. Körschen<sup>h</sup>, Florian Bergi<sup>i,j</sup>, Arild Folkvord<sup>i,j</sup>, U. Benjamin Kaupp<sup>h,k</sup>, Manfred Schartl<sup>l,m,n</sup>, and Leif Andersson<sup>a,o,p,2</sup>

<sup>a</sup>Science for Life Laboratory, Department of Medical Biochemistry and Microbiology, Uppsala University, 75123 Uppsala, Sweden; <sup>b</sup>Science for Life Laboratory, National Bioinformatics Infrastructure Sweden, Uppsala University, 75237 Uppsala, Sweden; <sup>c</sup>UR 1037 Fish Physiology and Genomics, INRAE, F-35000 Rennes, France; <sup>d</sup>State Key Laboratory of Developmental Biology of Freshwater Fish, College of Life Sciences, Hunan Normal University, Changsha 410081, China; <sup>e</sup>Fakultät für Biologie und Biotechnologie, Ruhr-Universität Bochum, 44801 Bochum, Germany; <sup>f</sup>Center for Marine and Molecular Biotechnology, Pilot National Laboratory for Marine Science and Technology (Qingdao), 266237 Qingdao, China; <sup>g</sup>College of Marine Life Sciences, Ocean University of China, 266100 Qingdao, China; <sup>h</sup>Molecular Sensory Systems, Center of Advanced European Studies and Research, 53175 Bonn, Germany; <sup>i</sup>Department of Biological Sciences, University of Bergen, 5020 Bergen, Norway; <sup>j</sup>Institute of Marine Research, 5817 Bergen, Norway; <sup>k</sup>Life & Medical Sciences Institute, University of Bonn, 53115 Bonn, Germany; <sup>l</sup>Developmental Biochemistry, Biocenter, University of Würzburg, 97074 Würzburg, Germany; <sup>m</sup>Hagler Institute for Advanced Study, Texas A&M University, College Station, TX 77843; <sup>n</sup>Xiphophorus Genetic Stock Center, Texas State University, San Marcos, TX 78666; <sup>o</sup>Department of Veterinary Integrative Biosciences, Texas A&M University, College Station, TX 77843; and <sup>p</sup>Department of Animal Breeding and Genetics, Swedish University of Agricultural Sciences, 75007 Uppsala, Sweden

Contributed by Leif Andersson, August 11, 2020 (sent for review May 20, 2020; reviewed by Henrik Kaessmann and Catherine L. Peichel)

The mechanisms underlying sex determination are astonishingly plastic. Particularly the triggers for the molecular machinery, which recalls either the male or female developmental program, are highly variable and have evolved independently and repeatedly. Fish show a huge variety of sex determination systems, including both genetic and environmental triggers. The advent of sex chromosomes is assumed to stabilize genetic sex determination. However, because sex chromosomes are notoriously cluttered with repetitive DNA and pseudogenes, the study of their evolution is hampered. Here we reconstruct the birth of a Y chromosome present in the Atlantic herring. The region is tiny (230 kb) and contains only three intact genes. The candidate male-determining gene *BMPR1BBY* encodes a truncated form of a BMP1B receptor, which originated by gene duplication and translocation and underwent rapid protein evolution. *BMPR1BBY* phosphorylates SMADs in the absence of ligand and thus has the potential to induce testis formation. The Y region also contains two genes encoding subunits of the sperm-specific Ca²⁺ channel CatSper required for male fertility. The herring Y chromosome conforms with a characteristic feature of many sex chromosomes, namely, suppressed recombination between a sex-determining factor and genes that are beneficial for the given sex. However, the herring Y differs from other sex chromosomes in that suppression of recombination is restricted to an ~500-kb region harboring the male-specific and sex-associated regions. As a consequence, any degeneration on the herring Y chromosome is restricted to those genes located in the small region affected by suppressed recombination.

sex determination | BMPR1 | CatSper | gene duplication | molecular evolution

The mechanisms underlying sex determination are astonishingly plastic. Particularly the triggers for the molecular machinery, which recalls either the male or female program, are highly variable and have evolved independently and repeatedly (1–8). Fish show a rich diversity of sex determination systems, including both genetic and environmental triggers (1, 9, 10). Adaptive hypotheses prevail in explaining the evolution of the various sex-determining triggers that guarantee a stable sex ratio of the population. Genetic sex determination is generally explained to be stabilized through the evolution of sex chromosomes, when sexually antagonistic selection links the genetic sex determinant to a locus that is beneficial to the same sex or even antagonistic to the opposite sex (11–15).

The Atlantic herring (*Clupea harengus*) is one of the most abundant vertebrates on Earth. Neither its karyotype (16) nor the male and female linkage map (17) indicated the presence of

sex chromosomes. Unexpectedly, our analysis of the Atlantic herring genome revealed a Y chromosome having a minute male-specific region (~230 kb), a small sex-associated region (~300 kb), and a large pseudoautosomal region (~30 Mb) corresponding to chromosome 8 in the reference assembly (17). Here we reconstruct the birth of this chromosome involving two gene duplications followed by rapid protein evolution of the candidate male-determining gene *BMPR1BBY*.

## Results

**Identification of the Sex-Determining Region.** The sex-associated region was identified by genome-wide association analysis using whole-genome sequencing data aligned to the reference genome (17). We compared SNP (single nucleotide polymorphism) allele

### Significance

Understanding the evolution of sex determination mechanisms and sex chromosomes is of fundamental importance in biology. Here we have reconstructed the evolution of the sex-determining region in the Atlantic herring. The region is small and contains only three intact genes. The candidate sex-determining factor *BMPR1BBY* is an evolutionary innovation in the herring lineage. It encodes a truncated form of a BMP type I receptor, which originated by gene duplication and underwent rapid protein evolution. The receptor has maintained its kinase activity and has the potential to induce development of testis. The other two genes in the sex-determining region, *CATSPERG* and *CATSPER3Y*, are male beneficial genes because they encode proteins predicted to be essential for sperm to fertilize the egg.

Author contributions: N.R., A.H., U.B.K., M.S., and L.A. designed research; N.R., J.C., A.H., M.E.P., F.H., C.F., O.W., C.-J.R., S.P., A.C., M.L., C.T., A.P., K.K., W.B., and H.G.K. performed research; F.B. and A.F. contributed new reagents/analytic tools; N.R., J.C., A.H., M.E.P., F.H., O.W., C.-J.R., M.L., C.T., A.P., K.K., W.B., H.G.K., U.B.K., M.S., and L.A. analyzed data; and N.R., J.C., A.H., M.E.P., U.B.K., M.S., and L.A. wrote the paper.

Reviewers: H.K., Heidelberg University; and C.L.P., University of Bern.

The authors declare no competing interest.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

¹N.R., J.C., A.H., and M.E.P. contributed equally to this work.

²To whom correspondence can be addressed. Email: leif.andersson@imbim.uu.se.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2009925117/-/DCSupplemental.
frequencies between eight males and eight females by $\chi^2$ analysis (Dataset S1); 765 SNPs spanning ~160 kb on chromosome 8 showed significant allele frequency differences ($P \le 1.0 \times 10^{-10}$) (SI Appendix, Fig. S1A). Analysis of SNP genotypes across the sex-associated region for the 16 fish revealed that males were heterozygous (XY), whereas females were homozygous (XX) for most sex-associated SNPs (SI Appendix, Fig. S1B). Thus, males are the heterogametic sex. An additional 17 Atlantic herring of unknown sex followed this pattern perfectly with individuals either being heterozygous or homozygous across the sex-associated region. The sequence of this region in the current assembly (17) represents the female haplotype (SI Appendix, Fig. S1B). In conclusion, herring chromosome 8 can now be split into chromosome (Chr) X and Y.

To characterize the male-specific region on ChrY, we examined the sequence of all unplaced scaffolds from the current

assembly (Methods) (17). Two unplaced scaffolds, 18 and 229, showed regions of high identity with the female-specific sequence on ChrX (SI Appendix, Fig. S1C and D). These regions were interrupted by male-specific sequences that are present in only one copy in males (Fig. 1A and B). We also analyzed data from six individually sequenced Pacific herring of unknown sex (18) and found that four had no coverage in the male-specific regions, whereas two had a single copy of the male-specific regions (Fig. 1B). The results provide evidence that the male-specific region evolved prior to the divergence of the Atlantic and Pacific herring that separated ~2 My before present (19).

We used the Xdrop method (20) for targeted resequencing of long DNA molecules (up to ~100 kb) from the borders between the male-specific and the flanking region. The alignment of these reads to our ChrY assembly confirmed a correct orientation and did not reveal any gaps in the alignment (SI Appendix, Fig. S1D).

<figure><img src="image_118.png" /><figcaption>**Fig. 1.** Assembly of chromosome Y in Atlantic herring. (A) Assembly and annotation of male-specific region (red), sex-associated region (orange), and flanking pseudoautosomal region (blue). Some genes (annotated as "_like") within the male-specific region show partial similarity with known proteins. LTR track, long terminal repeats. d, nucleotide sequence divergence between ChrX and Y in 500-bp windows. (B) Normalized copy number of sequenced individuals (eight female and eight male Atlantic herring and six Pacific herring of unknown sex) across the male-specific and sex-associated regions. Regions with high coverage marked with an asterisk overlap repetitive sequences. (C) Genetic differentiation ($F_{ST}$) in 50-kb overlapping windows between males (n = 8) and females (n = 8) across ChrY. The SDR is indicated by a vertical line. The male-specific region was not included in this analysis, since it is missing on ChrX. (Inset) $F_{ST}$ values in the near vicinity of the SDR.</figcaption></figure>
Our assembly of the ~230-kb male-specific region of ChrY shows that the male-specific sequence is split into two parts (Fig. 1A, in red) separated by an ~110-kb sex-associated region (Fig. 1B, in orange) that is well assembled and annotated (17). It contains 13 genes that are present on both X and Y, as indicated by equal sequence depth in males and females (Fig. 1B). The region showing strong genetic differentiation between males and females is remarkably sharp and only extends for ~300 kb from ~21.0 to ~21.3 Mb on ChrX (Fig. 1C), and the male-specific region together with the sex-associated region extends for ~500 kb on ChrY (SI Appendix, Fig. S1). These results imply that there is no meiotic recombination in the male-specific regions (red) and no or at least suppressed recombination in the sex-associated regions (orange) (Fig. 1A). The genes located in the region showing strong genetic differentiation between X and Y are listed in SI Appendix, Table S1. The nucleotide distance between X and Y in the sex-associated region located between the male-specific regions (d = 0.96%; Fig. 1A) is only 2.5-fold higher than the estimated average nucleotide distance between Atlantic and Pacific herring (~0.4%), implying either a divergence time of ~5 My or ongoing genetic exchange in the sex-associated region.

The male-specific region contains only three intact genes—*CATSPERG*, *CATSPER3Y*, and *BMPR1BBY*—that all are of high interest for male function (Fig. 1A). In addition, we identified other genes annotated ab initio, but these were not well supported by RNA-seq data from gonadal tissue; one of these, POL-like, is of retroviral origin. *CATSPERG* and *CATSPER3Y* encode subunits of the sperm-specific Ca²⁺ channel CatSper that is required for sperm function and male fertility in mouse, humans, and sea urchin (21–25). A phylogenetic tree analysis shows that the substitution rates of *CATSPERG* in Atlantic herring and its orthologs in other teleost are similar, consistent with a maintained function as a CatSper subunit (SI Appendix, Fig. S2). *CATSPER3* exists in two copies, *CATSPER3Y* in the male-specific region and a pseudoautosomal copy, *CATSPER3A*, ~1 Mb apart from the Y copy (Chr8, 22,259,380 to 22,268,493 bp). We estimate that *CATSPER3Y* arose by duplication about 15 Mya (SI Appendix, Fig. S3). Although the two *CATSPER3* genes in herring are closely related, *CATSPER3Y* is more well conserved

compared with paralogs in other species (SI Appendix, Fig. S3), suggesting that its function is maintained. By contrast, two findings suggest that *CATSPER3A* is nonfunctional: its exon–intron organization is not intact, and the dN/dS ratio of nonsynonymous (dN) to synonymous (dS) substitutions is close to 1.0, which is characteristic for pseudogenes.

The third gene in the male-specific region *BMPR1BBY* encodes a severely truncated form of BMPR1B (bone morphogenetic protein receptor 1B), a transmembrane protein that hosts an intracellular Ser/Thr kinase domain and is a member of the TGFβ receptor family. Components of TGFβ signaling control the development of male gonads (4). Sequence comparison of *BMPR1BBY* with the paralogs *BMPR1BA* and *BMPR1BB* on chromosomes 7 and 21, respectively, demonstrated that *BMPR1BBY* is more closely related to *BMPR1BB* (SI Appendix, Fig. S4). Transcriptome analysis of testis RNA revealed two isoforms with different 5' start sites (S- and L-BMPR1BBY) containing only coding sequences corresponding to exon 7 to 11 of the autosomal paralogs (SI Appendix, Fig. S4). In summary, *BMPR1BBY* is a truncated copy of BMPR1BB that has been translocated to chromosome 8; it is a strong candidate for being the male sex-determining factor.

We genotyped 200 individuals with known sex for the presence of *BMPR1BBY*. The high consistency between the X/Y genotype and sex phenotype (94%; SI Appendix, Table S2) demonstrates strong association to sex determination. Eleven out of 12 discrepancies between genotype and phenotype concerned XY females carrying *BMPR1BBY*, suggesting that the herring Y may not show full penetrance. Such incomplete penetrance of a young genetic sex determination system has also been observed in the northern pike (26).

**Tissue Expression of Genes from the Male-Specific Region.** We compared the tissue expression of *BMPR1BBY*, *CATSPER3Y*, *CATSPERG*, and the autosomal copies *BMPR1BB* and *CATSPER3A*, as well as two other TGFβ signaling genes, the anti-Müllerian hormone (*AMH*) and its cognate receptor (*AMHR2*), involved in vertebrate sexual development. No expression of *CATSPER3A* was detected. The comparison revealed that the

<figure><img src="image_119.png" /><figcaption>**Fig. 2.** Evolution of the ChrY-specific *BMPR1BBY* gene and protein. (A) Schematic representation of herring BMPR1BB and BMPR1BBY proteins showing the architecture of BMP receptors and highlighting the absence of the extracellular, transmembrane and GS-box domains in the truncated BMPR1BBY. (B) dS and dN phylogenetic trees of *BMPR1BB* (Herring A) and *BMPR1BBY* (Herring Y) compared to *BMPR1BB* in other species. (C) Crystal structure of the intracellular domain of human BMPR1B (PDB 3MDY, amino acids 206 to 500, red) and a model of herring BMPR1BBY (blue). Yellow color represents the N-terminal extension of L-BMPR1BBY compared to S-BMPR1BBY. Orange color indicates amino acid substitutions between BMPR1BB and BMPR1BBY. N, N terminus; C, C terminus.</figcaption></figure>
three Y-specific genes were only expressed in males as expected (SI Appendix, Fig. S5). *BMPR1BBY* and, surprisingly, also *CATSPER3Y* were expressed in both testis and brain, whereas *CATSPERG* showed testis-specific expression. The autosomal copy of *BMPR1BB* was more broadly expressed with the highest level in gonads; female gonads showed higher transcript levels than male gonads (SI Appendix, Fig. S5B). As expected, *AMH* and *AMHR2* were highly expressed in male and female gonads (SI Appendix, Fig. S5 D and F), consistent with their important roles in sex determination and gonad differentiation.

**Rapid Protein Evolution of BMPR1BBY.** The predicted structure of the autosomal herring BMPR1BB and the canonical vertebrate BMPR1 receptor is similar (Fig. 2A) (27). The *BMPR1BBY* transcripts are lacking exons 1 to 6 of the autosomal *BMPR1BB* (SI Appendix, Fig. S4) that encode the extracellular, transmembrane, and part of the intracellular domain (Fig. 2A). These domains mediate ligand binding, receptor activation, and membrane localization and entail the phosphorylation sites of the GS loop. Amino acid alignment of BMPR1BBY and BMPR1BB from herring and other fish reveals more than 40 changes that have occurred since BMPR1BBY arose (SI Appendix, Fig. S4C). The rate of dS and dN substitutions between BMPR1BBY and BMPR1BB is 0.12 and 0.07, respectively, yielding a dN/dS ratio of 0.58, a dramatic and highly significant (P < 10⁻²¹) increase compared with dN/dS ratios of 0.01 to 0.02 in other branches of the tree (Fig. 2B and SI Appendix, Fig. S6). We estimate that BMPR1BBY arose by duplication about 20 Mya (SI Appendix, Fig. S6).

We built a structural model of the intracellular domain of herring BMPR1BBY and BMPR1BB based on human BMPR1B (Fig. 2C). Both herring BMPR1BB protein folds are essentially identical to that of the human paralog (Cα root-mean-square deviation of 0.008, 0.077, and 0.073 Å for BMPR1BB, L-BMPR1BBY, and S-BMPR1BBY, respectively; SI Appendix, Fig. S7 and Table S5), suggesting that the truncated form may still display Ser/Thr kinase activity.

To conclude, *BMPR1BBY* shows an accelerated protein evolution typical for sex-determination genes (28) but appears to maintain the BMPR1B protein fold.

**Physiological Activity of BMPR1BBY.** TGFβ signaling by the anti-Müllerian hormone (AMH) system is critically involved in sex determination (29). BMPR1BB, a type-I TGFβ receptor, is activated by ligand-bound type-II receptor and, in turn, phosphorylates the transcription factor SMAD. AMH-induced signaling occurs primarily via SMAD5 phosphorylation and to some extent via SMAD1 and SMAD8 (30). To probe the kinase activity of truncated BMPRBBY, we established an in vitro reporter system for SMAD activation (Fig. 3A and Methods). The system consists of a Gal4 transcriptional activator coupled to the transactivation domains of SMAD1 and SMAD5, which are specific for AMH signaling (31), and a UAS-Luc reporter driving luciferase expression. We tested three different scenarios: 1) the canonical AMH-signaling machinery including herring AMH, AMHR2, and BMPR1BB; 2) the canonical AMH machinery for which BMPR1BB has been replaced by either S-BMPR1BBY or L-BMPR1BBY; and 3) the canonical AMH machinery supplemented by either S-BMPR1BBY or L-BMPR1BBY (Fig. 3 B–E). Phosphorylation levels of SMAD1 did not differ significantly when BMPR1BB was replaced by either S-BMPR1BBY or L-BMPR1BBY (Fig. 3 B and D). By contrast, S-BMPR1BBY is fourfold more active than BMPR1BB in SMAD5 phosphorylation (Fig. 3C). Both SMAD1 and SMAD5 activity were enhanced when both BMPR1BB and S-BMPR1BBY were present (Fig. 3 B and C). No significant increase of SMAD1 and only low increase of SMAD5 phosphorylation were recorded after addition of L-BMPR1BBY (Fig. 3 D and E). Thus, S- and L-BMPR1BBY can substitute for BMPR1BB in activating SMAD1 and SMAD5.

Next, we examined whether S- and L-BMPR1BBY synergize with canonical AMH signaling. To this end, AMHR2, BMPR1BB, and S- or L-BMPR1BBY were tested either in the presence or absence of the AMH ligand. While canonical AMH signaling (AMH, AMHR2, and BMPR1BB) was clearly stimulated by AMH (SI Appendix, Fig. S8A), AMH did not modulate S- or L-BMPR1BBY-mediated signaling when associated with AMHR2 (SI Appendix, Fig. S8 B and C). Moreover, S- or L-BMPR1BBY-mediated signaling via SMAD1 phosphorylation was unaffected when AMH was expressed without its cognate receptor AMHR2 (SI Appendix, Fig. S8 D and E), indicating that BMPR1BBY-mediated signaling is not further modulated by the AMH ligand. Finally, in the absence of AMH, the presence or absence of AMHR2 did not affect S- or L-BMPR1BBY-mediated signaling (Fig. 3F), while differential synergistic effects were observed for BMPR1BBYs when expressed together with the whole AMH machinery (AMH ligand, AMHR2, and BMPR1BB; Fig. 3 B–E).

In summary, BMPR1BBY has maintained its enzymatic activity despite the severe truncation and the accelerated evolution at the protein level. S- and L-BMPR1BBY-mediated signaling requires neither the presence of the AMH ligand nor the AMHR2 receptor or the autosomal BMPR1BB, yet is able to enhance canonical AMH signaling. Thus, the Y-linked receptor has the enzymatic features to act as an upstream regulator of male gonad development.

**CatSper3Y and CatSperG Are Essential Components of the Herring Sperm Proteome.** The location of two CatSper genes on the herring Y chromosome is intriguing because CatSper has been thought to be lost in several metazoan lineages including fish (32). We revisited this issue and identify here CatSper genes in many but not all fish genomes (Fig. 4 A and B and Dataset S2). CatSper is a uniquely complex Ca²⁺ channel that consists of at least nine different subunits CatSper1–CatSper4, CatSperB, CatSperG, CatSperD, CatSperE, and CatSperZ. The Atlantic herring genome contains all CatSper genes except CatSperZ (33) (Dataset S2). Using mass spectrometry (MS), seven CatSper subunits including CatSper3 and CatSperG were identified in mature herring sperm (Datasets S3 and S4). Signaling pathways control CatSper and, thereby, sperm motility, navigation, and fertilization from marine invertebrates to mammals (34–36). The CatSper channel is gated open by membrane voltage and alkaline pH. Three signaling events are key to CatSper activation in sperm of the sea urchin *Arbacia punctulata*: hyperpolarization of the membrane potential Vₘ by a K⁺ channel CNGK (37, 38) followed by repolarization via hyperpolarization-activated and cyclic nucleotide-gated (HCN) channels, so-called pacemaker channels (39, 40), and alkalization by a sperm-specific voltage-gated sodium/proton exchanger (sNHE or SLC9C1) (41). The activity of the SLC9C1 exchanger and HCN channel is modulated by cAMP, which is synthesized by a soluble adenylate cyclase (sAC or *acyl*10 gene). Remarkably, all these signaling proteins upstream of CatSper are present in herring sperm, suggesting that CatSper and Ca²⁺ also control motility of herring sperm (Fig. 4C and Dataset S4). Two functionally important variations are striking. First, herring sperm carry an HCN-like (HCNL) channel that is distinct from classic HCN channels and falls into a new gene subfamily (42). HCNL1 in zebrafish sperm, unlike classic HCN channels, is proton-selective and not cAMP-sensitive (42). Although herring sperm carry the HCNL2 isoform that has not yet been characterized functionally, structural features underlying proton permeation are conserved between HCNL1 and HCNL2 (34). Second, the zebrafish CNGK channel is also gated not by cyclic nucleotides but by alkaline pH, indicating that the cellular messengers in fish sperm are protons rather than cAMP or cGMP (42, 43). These results suggest that herring sperm employ a signaling pathway for fertilization that combines signaling proteins and mechanisms from marine invertebrates, zebrafish, and mammals. As more signaling molecules
<figure><img src="image_120.png" /><figcaption>**Fig. 3.** Analysis of BMPR1BBY signaling. (A) Schematic of the UAS-Gal4-Smad-AD reporter assay for quantifying the degree of SMAD1 and SMAD5 phosphorylation upon application of different combinations of TGFβ ligands and receptors (see Methods for detailed information). In absence of AMHR2 or BMPR1BBY-induced stimulations, the Gal4-Smad-AD fusion protein remains in the cytoplasm and therefore does not activate transcription of the UAS-4 promoter driving expression of the luciferase. In this situation, only background luciferase expression levels are recorded and set a blank value. In presence of an active type I receptor, the Gal4-Smad-AD fusion proteins are selectively phosphorylated by the receptor through the Smad-AD domain and then translocate to the nucleus. Here they bind and transactivate the UAS-4 promoter through its Gal4 DNA binding domain, which leads to transcription of the luciferase gene. The resulting luciferase activity thus strictly correlates to the degree of activation (phosphorylation) of the Smads, which in turn is a measure of the specific upstream signals elicited by the various combinations of ligands and receptors. AD, specific transactivation domains of the different SMADS. (B–E) Normalized luciferase intensities relative to controls after transfection of different combinations of herring AMH, AMHR2, BMPR1BB, and S- and L-BMPR1BBY. (F) Normalized luciferase intensities relative to controls using the SMAD1 reporter after transfection of different combinations of herring AMH, AMHR2, BMPR1BB, and S- and L-BMPR1BBY.</figcaption></figure>

are identified, the herring data strengthen an emerging picture of evolutionary significance: sperm preferentially adopt unique, sperm- and species-specific signaling molecules and pathways that are tuned to their respective habitat of reproduction (35).

In summary, we conclude that the herring Y-specific region contains genes encoding two proteins, CatSperG and CatSper3Y, expected to be essential for sperm function based on the functional conservation of the CatSper channel in both invertebrates and vertebrates.

Discussion

BMPR1BBY is a very strong candidate sex-determining factor in the Atlantic herring. It is clearly the best candidate out of the only three genes in the male-specific region because of the well-established role for TGF/BMPR signaling in sexual development (29). Its presence is strongly correlated with the male phenotype, it shows an accelerated protein evolution, and it can replace the AMH/BMPR1B/AMHR pathway for SMAD phosphorylation.
<figure><img src="image_121.png" /><figcaption>**Fig. 4.** Characterization of CatSper genes and the Ca²⁺ signaling pathway in sperm. (A and B) Multiple independent gene losses of CatSper1 to CatSper4 and CatSperB, CatSperG, CatSperD, and CatSperE in fish. Cladogram of major phylogenetic groups based on Betancur et al. (70). Colored branches indicate the presence of at least three of each the CatSper1 to CatSper4 and CatSperB, CatSperG, CatSperD, and CatSperE in all (green for both and blue for CatSper1 to CatSper4 only), none (orange), or some (green and orange for both, blue and orange for CatSper1 to CatSper4 only) species of the clade. Number of species with CatSper genes/all species investigated. Typical members of the respective clade are shown on the right. (A) Full tree including all fish species with searchable genomes. Clade groupings are simplified: early diverging bony fish comprise several groups—*Cladistia*, *Chondrostei*, coelacanths, eels, and *Holostei*. Note the near-complete retention of CatSper in cartilaginous fish and salmon-related fish (*Protacanthopterygii* and *Stomiatii*). (B) Presence and absence of CatSper in the *Otomorpha* clade. Note the almost complete retention of CatSper1 to CatSper4 in *Cyprinidae*, whereas all three sister clades have lost CatSper1 to CatSper4. Within *Clupeiformes*, four of five species have retained CatSper. (C) Comparison of signaling pathways in sperm of the sea urchin (*Arbacia punctulata*), Atlantic herring (*C. harengus*), and zebrafish (*D. rerio*). GC, chemoreceptor guanylate cyclase; CNGK, K⁺-selective cyclic nucleotide-gated channel; sNHE, sperm specific Na⁺/H⁺ exchanger SLC9C1; HCN1 and HCN2, hyperpolarization activated and cyclic nucleotide-gated channels; CatSper, sperm-specific Ca²⁺ channel; NCKX, Na⁺/Ca²⁺/K⁺ exchanger; PMCA, plasma membrane Ca²⁺-ATPase; HCNL1, HCN-like 1 channel; HCNL2, HCN-like 2 channel.</figcaption></figure>

The dN/dS ratio of 0.58 for *BMPR1BBY* is dramatically higher than the dN/dS = 0.01 to 0.02 for the autosomal paralogs (*SI Appendix*, Fig. S6). However, a dN/dS ratio of 0.58, well below 1.0, may be explained by relaxed purifying selection or a combination of positive selection at some sites and purifying selection at other sites. Our findings that the structural modeling predicts a conserved protein structure as well as its conserved kinase activity strongly suggest that the latter is the correct explanation. An important topic for future studies will be to further characterize the expression pattern of *BMPR1BBY* and in particular explore if it is expressed during early testis development at the sex determination stage. However, such experiments are challenging in a pelagic fish that does not reproduce in captivity.

The herring Y chromosome provides a unique opportunity to reconstruct the birth of a sex chromosome due to the relatively small size of the male-specific region that is not overloaded with

repetitive DNA and other junk material like most other sex chromosomes. The evolutionary process must have started with the emergence of a truncated copy of *BMPR1BB* that was translocated to chromosome 8 (Fig. 5). The truncation assigned a new protein function to the gene. At the same time, because of the loss of its 5' end including the promoter and intron 1, *BMPR1BBY* must have lost the expression regulation of the autosomal gene and most likely acquired a new spatial and temporal expression pattern required to trigger male gonad development. Such transcriptional rewiring of emerging sex-determining genes, in particular if derived from gene duplication of a downstream gonadal factor, has been postulated to be a hallmark of sex chromosome evolution (4, 44).

The insertion of the duplicated fragment from chromosome 21 into chromosome 8 created a region of nonhomology that cannot recombine in meiosis. Thus, a proto Y was formed and maintained
by suppressed recombination ab initio. It is likely that *CATSPER3* and *CATSPERG* were already colocalized on the ancestral chromosome 8 when *BMPR1BBY* arrived on the scene because they are located on the same chromosome in many teleosts (*SI Appendix*, Table S3). In goldfish, the two genes are <1 Mb apart. Possibly the next step was duplication of *CATSPER3* and its translocation near *BMPR1BBY* during or after the duplication event (Fig. 5). The synteny in that region is highly conserved across vertebrates, including *CATSPER3* and the homeobox gene *PITX1*, although CatSper genes have been lost in birds and many teleosts (*SI Appendix*, Fig. S9, and Fig. 4). Notably, this region is more rearranged in the herring than in other vertebrates (*SI Appendix*, Fig. S9), suggesting that chromosomal rearrangements contributed to the evolution of the Y-specific region. A third event was the incorporation of *CATSPERG* into the Y-specific region and its loss from a pseudoautosomal location. Finally, *CATSPER3A* underwent the fate of most duplicated genes, namely, degeneration, because it was redundant and *CATSPER3Y* benefited from its tighter linkage to the male sex-determining gene. By this process, a minimal sex-determining region (SDR) has evolved in the Atlantic herring; the region does not recombine with chromosome X and contains only three functional genes, a candidate male-determining factor and two genes essential for sperm function.

It is possible that the SDR in herring as present today is much younger than the emergence of the BMPR1BBY gene about 20 Mya (SI Appendix, Fig. S6), as indicated by the modest sequence divergence in the sex-associated region which rather suggests a divergence time of a few million years (Fig. 1A). However, it is possible that there is ongoing genetic exchange (recombination or gene conversion) between X and Y in this region. In fact, one of the Pacific herring chromosomes included in Fig. 1B may represent a recombinant haplotype because it appears X-like over the *BMPR1BBY-CATSPER3Y* region but Y-like over *CATSPERG*. It will be of considerable interest to further explore the evolution of this SDR region in a closely related clupeid like the sprat (*Sprattus sprattus*).

The association of *BMPR1BBY* with *CATSPER3Y* and *CATSPERG* in the herring SDR is consistent with one of the most characteristic features of sex chromosomes, namely, suppressed recombination between a sex-determining factor and genes beneficial for that sex (11–15). However, the Atlantic herring Y differs from most previously described sex chromosomes as this suppressed recombination has not spread along the chromosome. The region showing strong genetic differentiation between males and females is remarkably small and extends only over ~300 kb despite the existence of the herring Y for at least a few

million years. As a comparison, the Y chromosome in the threespine stickleback, which is less than 26 My old, has a non-recombinating region of 17.5 Mb that shows extensive degeneration (45). Unlike the stickleback Y chromosome, the herring Y does not lack any gene present on the X chromosome; it only carries three additional genes that are likely essential for male development and function. Any degeneration that may have occurred is expected to be restricted to the genes located in the sex-associated region showing sequence divergence between X and Y (SI Appendix, Table S1).

## Methods

**Whole-Genome Sequence Analyses.** We analyzed whole-genome sequence (WGS) data of 16 individuals with known sex, collected from Atlantic and Baltic with on average ~12X sequence coverage (19). The reads were aligned to two assemblies of the Atlantic herring genome (17, 19), and high confidence SNPs were called as described (17, 19). A genome-wide association analysis for the detection of sex-associated regions was done by comparing allele frequency between males and females with a $\chi^2$ test using read counts extracted from sites called in the older assembly (19). The depth of coverage across the chromosome-based genome assembly (17) was extracted using GATK:DepthOfCoverage (version 3.5.0) (46). The data were converted to normalized depth 1-kb windows, and we compared the depth between males and females. No Animal Care and Use Committee approval was required because the herring were not harvested specifically for research purposes.

**Assembly of the Chromosome Y Region.** The assembly of the Y region is based on two unassigned scaffolds (18 and 229) and parts from an earlier version of our current PacBio-assembly of the herring genome (17), which is based on a male Atlantic herring, but only the ChrX region was included in the previously reported assembly (17). The Y assembly presented here was validated by investigating uniformity of depth of coverage of individual PacBio reads across the assembled sequence.

We used eight males and eight females to estimate the genetic differentiation between males and females by Weir–Cockerham $F_{ST}$ (47) using vcftools v0.1.16 (48) in 50-kb overlapping windows across chromosome Y.

**Annotation of the Male-Specific Region.** We identified male-unique regions on unassigned scaffolds from the current PacBio-based herring assembly which was not annotated by the Ensembl annotation pipeline (https://www.ensembl.org) because the male-unique region was not included in this assembly (17). Thus, we annotated these regions by the annotation pipeline implemented by National Bioinformatics Infrastructure Sweden based on the Maker package (version 3.01.02) (49) (https://github.com/NBISweden/GAAS, version 2019-07-11). We first generated an evidence-built annotation based on protein similarity using Uniprot (50) (downloaded August 11, 2016). We also generated an ab initio built annotation by first training Augustus (version 2.7) (51) using 1,213 genes annotated by Ensembl (17). Then we combined the two annotations to generate the final annotation. We manually curated *BMPR1BBY* by comparing it with the autosomal copy on chromosome 21 and generated the final gene structure. In addition, we

<figure><img src="image_122.png" /><figcaption>**Fig. 5.** Reconstruction of the evolution of sex chromosomes in Atlantic herring. The evolution of ChrY occurred in four steps, but the exact order of steps 2 to 4 is not yet known. 1) Duplication and translocation of *BMPR1BB* from Chr21 to Chr8. 2) Duplication and translocation of *CATSPER3* within Chr8/ChrY. 3) Incorporation of *CATSPERG* in ChrY and loss from ChrX. 4) *CATSPER3A* becomes pseudogenized or evolved a new function. MYA, million years ago.</figcaption></figure>
used GenScan (52) to do another round of ab initio annotation which was
merged with the Maker annotation.

To validate the structure of the gene model we conducted 3' RACE and Nanopore whole-transcriptome and gene-specific cDNA sequencing. Total RNA was prepared from the testis of an adult Atlantic herring using the RNeasy Mini Kit (Qiagen). The 3' RACE reaction was performed with 1 µg RNA using the FirstChoiceTM RLM-RACE Kit (Thermo Fisher Scientific). Nested RACE PCRs were performed in a 25-µL reaction containing KAPA HiFi Fidelity Buffer with MgCl₂, 0.32 mM dNTPs, 0.4 µM each of the forward and reverse primer, 0.8 U KAPA HiFi DNA Polymerase (Kapa Biosystems) and 1 µL of the cDNA or Outer RACE PCR product (1:5 dilution) as PCR template. Amplification was carried out following the program: 95 °C for 5 min, 35 cycles of 98 °C for 20 s, 58 °C for 30 s and 72 °C for 1 min, and a final extension of 10 min at 72 °C.

**Targeted Resequencing Using Xdrop.** To confirm the orientation of the male-specific region in relation to the flanking pseudoautosomal region we used the Xdrop method for targeted sequencing of long DNA molecules (20). Primer pairs for the amplification of fragments from the male-specific region in the near vicinity of the border to the pseudoautosomal region were designed. The following primer pairs were used to amplify a fragment in the vicinity of *CATSPER3Y*: for dPCR, Fw-TGCACTAGTGGGTGACTACCAA and R-TGTAAGACACCAGAATTTCCCCAC (product size 173 bp, Tm = 76.25 ± 0.08 °C), and for qPCR, Fw-TACTCAGATGGCCGACGAAATCTTTA and R-GCGTCAAAGACGTTATATCCGC (product size 111 bp, Tm = 77.88 ± 0.11 °C). For the amplicon from the scaffold18 region, the following primers were used: for dPCR, Fw-ACCGGCTTTATCCATCCGTC and R-AGTGACGAAAACGAAACC GTC (product length 165 bp, Tm = 76.97 ± 0.08 °C), and for qPCR, Fw-CTTCGGATCTGCACGATTCAC and R-TCAGAAGACCGAGACTGGA (product length 127, Tm = 80.98 ± 0.11 °C). The reactions were run in a final volume of 10 µL, of which 5 µL were of dPCR mix, 0.5 µL of qPCR dye, 0.4 µL of 10 µM primer pair, and 0.1 µL of water. The enrichment and Nanopore sequencing of the genomic regions was performed by Samplix services (https://samplix.com/).

**Genotyping Sex-Specific Markers.** We genotyped 200 individuals for the presence/absence of *BMPR1BBY*. These samples had accurate sex records and were collected from two localities, Björköfjärden, Kattegat and Brofjorden, Skagerrak (19). The following primer pair was designed using primer3 with default settings, bby-fw CCAGACTGACTTGTGGTTCC and bby-R CTTCCCAGC AACAAGATGG, giving a 132-bp product with a melting temperature of 85.89 ± 0.07 °C (average and SD, respectively). The PCR reactions were composed of 5 µL SYBR Green PCR Master Mix (Applied Biosystems by Thermo Fisher Scientific), 1 µL of primer pair working solution at a concentration of 4 µM, and 4 µL of the DNA samples at a final concentration of 1 ng per reaction. Control reactions with no template were included at all points to spot any primer dimer formation. PCR was done using a QuantStudio 6 Flex instrument (Applied Biosystems by Thermo Fisher Scientific) with the thermal profile 2 min at 50 °C and 2 min at 95 °C, followed by 40 cycles of 15 s at 95 °C, 30 s at 57 °C, and 30 s at 72 °C. Melting curves were generated with incubation for 15 s at 95 °C, 1 min at 60 °C, and 15 s at 95 °C. The ramp rate between the two last steps was 0.05 °C/s. Primer sequences are also in SI Appendix, Table S4.

**Nanopore Sequencing.** Total RNA (500 ng) from testis was converted into cDNA following the Oxford Nanopore Technologies (ONT) protocol for rapid cDNA sequencing with the SuperScript IV reverse transcriptase (Invitrogen). To improve yield of full-length cDNA the incubation times were extended to 50 min at 50 °C and 20 min at 42 °C. The cDNA was amplified with LongAmp polymerase (NEB) and click primers (cPRM, ONT) for 14 cycles with 6 min elongation. Sequencing was done as partial runs on R9.4 and R9.5 flow cells as 1D reads. *BMPR1BBY* and *CATSPERG* were amplified from testis cDNA with LongAmp and internal primers and in combination with the ONT forward and reverse adaptor primers in order to capture the full transcript. PCR products were sequenced using the LSK309 1d² ligation sequencing protocol (ONT).

**Genomic Alignments.** Similarity between *CATSPER3Y* and *BMPR1BBY* and their presumed homologous autosomal loci was assessed by performing local alignments of the genomic nucleotide sequence, including flanking regions, of the two Y-specific gene copies against genomic sequence fragments, extracted from assembly version Ch_v2.0.2 (GCA_900700415.1) that covered the potentially orthologous autosomal genes. The targets were identified from the Ensembl-provided annotation of the herring genome (Database version 99.202). The alignments were performed using Clustal

Omega version 1.2.3 (53) using default parameters and analyzed in R (54) using the package "ape" version 5.3 (55).

We compared the gene content in the sex-associated region on X and Y using Clustal Omega version 1.2.3 (53) implemented in "ape" R package version 5.3 (55) as described above.

**Phylogenetic Analysis by Maximum Likelihood (PAML).** First, we collected a reference set of *CATSPER3* and *BMPR1BB* genes from the ENA database. After manual inspection of the possible transcripts, we retained the option for each species that most closely resembled the herring Y copy. Second, in-frame nucleotide alignments spanning the extent of the Y-specific copies were generated for both genes. For *BMPR1BB*, which is highly conserved, the nucleotide alignments were already in frame, as verified by the "trans" function from "ape" and subsequent comparison with the reference sequences. For *CATSPER3*, which is more variable, PAL2NAL (56) was used to translate the amino acid alignment to a nucleotide equivalent using the reference cDNA sequences. In both cases, ambiguous positions were discarded during the rate analysis. Next, we extracted the taxa represented in each alignment from the Fish Tree of Life (57) to generate unrooted guide trees for the PAML analysis. These trees were edited to split the herring entry into a Y-specific and autosomal copy, and for *CATSPER3*, the coelacanth was added as an additional unresolved branch. Finally, the alignments and guide trees were fed into the "codeml" module of PAML version 4.9 (58) and run under three different models: 1) All branches evolve at the same rate (NULL model). 2) Only the "herring Y-" (for *BMPR1BB*) or the "herring A-" (for *CATSPER3*) branch has a unique rate (two-rates model). 3) Both the "herring Y" and "herring A" branches have a unique rate (three-rates model).

The models were then compared using a likelihood ratio test statistic, which was evaluated against a $\chi^2$-distribution with degrees of freedom equal to the difference in estimated parameters between the compared models. The PAML control files, guide trees, and alignments used are found at https://github.com/nimarafati/Atlantic_herring_SDR.

We used *dS* values from the PAML analysis to estimate divergence times by extracting fossil-calibrated LCA estimates from the Fish Tree of Life (57) and dividing the time interval proportionally to the lengths of the terminal and internal branches.

**Structure Modeling.** Structural prediction of the intracellular domain of BMPR1BB, full-length L-BMPR1BBY, and full-length S-BMPR1BBY was performed using the RaptorX server (59), which gave the best models using Protein Data Bank 3MDY (human BMPR1B) (60) as the template. Structure alignments and protein representations were generated using the PyMOL Molecular Graphics System (Version 1.3, Schrödinger, LLC).

**Preparation and Mass Spectrometric Analysis of Flagella from Herring Sperm.** The preparation of flagella from *C. harengus* sperm was as described (39) using Herring Ringer's (HR) solution instead of artificial seawater. Briefly, sperm were washed two times with HR (in mM: 206 NaCl, 7.2 KCl, 2.1 CaCl₂, 3.1 MgCl₂, 10 Hepes/NaOH, pH 7.6) (61). Hepes was used for buffering instead of NaHCO₃. Sperm suspension was sheared 20 times by centrifugation for 30 s at 75 × g and 4 °C through a 40-µm mesh of a cell strainer (BD Biosciences). Flagella were washed in HR and stored as pellet.

Flagella protein (20 µg) was subjected to GeLCMS as described (62) with the following modifications: excision of 10 gel slices, use of Proteome Discoverer Rev. 2.2, *C. harengus* protein sequences from National Center for Biotechnology (NCBI) release GCF_900700415.1 (date April 16, 2019) supplemented with additional sequences (Dataset S4), recalibration of precursor MS spectra, precursor mass error tolerance 8 ppm, MS/MS mass error tolerance 0.6 Da, oxidation of Met, phosphorylation at Ser, Thr, Tyr, and N-terminal pyro-Glu as variable, Cys-Propionamide as static modification.

**Bioinformatic Analysis of CatSper Genes.** All fish genomes available from NCBI (https://www.ncbi.nlm.nih.gov/genome/) as of March 11, 2020, were analyzed for the presence of highly conserved β actin as positive control. The whole-genome shotgun databases (WGS) for these genomes were searched with TBlastN (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=tblastn) for β actin, using the *Danio rerio* (zebrafish) sequence as query. Genomes which were unsearchable or yielded poor hits were excluded. The research resulted in 271 genomes, which were evaluated for the presence of *CatSper1* to *CatSper4* and *CatSperB*, *CatSperG*, *CatSperD*, and *CatSperE*. We then performed TBlastN searches in the WGS databases of fish genomes using previously annotated protein sequences for *CatSper1* to *CatSper4* and *CatSperB*, *CatSperG*, *CatSperD*, and *CatSperE* of *Lepisosteus occulus* (spotted gar), *Latimeria chalumnae* (coelacanth), *Callorhinchus milii* (elephant shark), *Salmo salar* (Atlantic salmon), *Esox lucius* (northern pike), and *Carassius*
auratus (goldfish). The search was recursive, whenever more closely related sequences were identified. Candidates were evaluated according to e value until outgroup genes were found; a rigid cutoff by e value was ineffective. Candidate genes were predicted by genewise (https://www.ebi.ac.uk/Tools/psa/genewise/) (63) using closely related validated sequences as template. To validate candidates, we employed phylogenetic tree analysis. Trees were constructed from the candidate genes, the query sequences and other validated CatSper genes (e.g., *Mus musculus* [mouse], for full list, see Dataset S2). For CatSper1 to CatSper4, representative members from related ion channel families were used as outgroup (CACNA and SCN channels of mouse and zebrafish, TPCN and TRP channels of zebrafish). CatSperB, CatSperG, CatSperD, and CatSperE served as each other's respective outgroup. Sequences were aligned using the MAFFT (multiple alignment using fast Fourier transform) tool version 7 (https://mafft.crcb.jp/alignment/server/) (64, 65) with the iterative refinement strategy called E-INS-I strategy and otherwise default parameters. Gap Strip Squeeze (https://hcv.lanl.gov/content/sequence/GAPSTREEZE/gap.html) was used to remove regions of the multiple sequence alignments which contained gaps in more than 90% of the sequences. Trees were subsequently calculated using the PhyML online tool, which is a phylogenetic tool that employs a maximum likelihood (ML) algorithm (http://www.atgc-montpellier.fr/phylm/) (66). Tree improvement was set to SPR (subtree pruning and regrafting), and branch support was set to $χ^2$-based aLRT (approximate likelihood ratio test). Otherwise the default settings were unchanged. Branch support for single CatSper1 to CatSper4 clades was always maximal; thus, CatSper1 to CatSper4 candidates that fell into these clades were considered validated. Branch support for outgroup clades (CACNA, SCN, and TRP) was maximal, and candidates which were sorted to the outgroup genes were discarded. For CatSperB, CatSperG, CatSperD, and CatSperE, branch support was always close to maximal; thus, candidates falling into these groups were considered validated.

**RNA Preparation and Expression Analysis.** Total RNA was extracted from the gonad, heart, spleen, kidney, gills, intestine, hypothalamus and sacillus vasculosus (BSH), and brain without BSH (brain) of three adult male and three female Atlantic herring using the RNeasy Mini Kit (Qiagen). RNA was reverse transcribed into cDNA with a High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific). Real-Time PCR with SYBR Green chemistry was performed for **BMPR1BB**, **CATSPER3Y**, **CATSPERG**, **AMH**, and **AMHR2** in 10 µL reactions of SYBR Green PCR Master Mix (Thermo Fisher Scientific), 0.3 µM primers, and cDNA template, with a program of an initial denaturation of 10 min at 95 °C followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. TaqMan Gene Expression assay (Thermo Fisher Scientific) containing 0.30 µM primers and 0.25 µM **BMPR1BBY** TaqMan probe (Integrated DNA Technologies) was performed to compare the relative expression levels of **BMPR1BBY**. **ACTIN** was used as the housekeeping gene in all expression assays. All of the primers and probes are listed in *SI Appendix, Table S4*.

The expression of **CATSPER3A** in gonads and brain was investigated by RT-PCR. Briefly, cDNA was produced with the High Capacity cDNA Synthesis kit (Thermo Fisher Scientific), without RNase inhibitor. As this kit does not contain any oligo-dT primer, that was added to the reaction master mix to a final concentration of 5 mM [Oligo(dT)₁₈ primer; Thermo Fisher Scientific]. Primer pairs were designed with the NCBI tool (**SI Appendix, Table S4**). RT-PCRs were run with the AmpliTaq Gold DNA Polymerase with Buffer II and MgCl₂ (Thermo Fisher Scientific) together with dNTPs (dNTP Mix, 5 mM each) at a final concentration of 200 mM each. Primer pairs were added to a final concentration of 0.4 mM each, and the working concentration of DNA polymerase was 2 U/reaction. The total sample amount was 40 ng in a final

volume of 10 µL. The reactions were incubated according to the following profile: 95 °C for 10 min, 35 cycles of 95 °C for 15 s, 57 °C for 30 s, 72 °C for 1 min, followed by 72 °C for 5 min. All reactions were checked by gel electrophoresis. No detectable amplification was observed for either primer pair.

**Transfection Experiments.** Full-length clones for **AMH**, **AMHR2**, and autosomeal **BMPR1BB** as well as **L-BMPR1BBY** and **S-BMPR1BBY** clones were designed including a fish consensus Kozak sequence (67), synthesized and cloned (GenScript) into the expression vector pcDNA3.1 (-) (Thermo Fisher Scientific).

Medaka spermatogonial stem cells (Sg3) were cultured as previously described (68). For transfection, cells were grown to 80% confluency in six-well plates and subsequently transfected with 5 µg expression vectors using the FuGene (Roche) reagent as described by the manufacturer.

For monitoring differential activation of Smads upon selective AMH, **AMHR2**, **BMPR1B**, or S- or L-BMPR1BBY expression, Sg3 cells were seeded in six-well plates and cotransfected with a combination of four kinds of plasmids (Fig. 3A): 1) an expression plasmid which encodes a fusion protein (31, 69) (Smad1-AD-Gal4-DBD or Smad5-AD-Gal4-DBD; 300 ng per well) that will translocate to the nucleus upon phosphorylation by the different TGFs and in return transactivate the UAS-4 promoter through its Gal-4 DNA binding domain; 2) a reporter plasmid, which codes for luciferase under the control of a minimal promoter, which contains UAS sequences (UAS-luc, *Firefly* luciferase; 300 ng per well); 3) plasmids coding for the different ligands and receptors to be tested for signaling activity [either pcDNA3.1(-)-AMH, pcDNA3.1(-)-AMHR2, pcDNA3.1(-)-BMPR1BB, pcDNA3.1(-)-S-BMPR1BY, pcDNA3.1(-)-L-BMPR1B or a control plasmid (empty pcDNA3.1(-)); 400 ng per well]; and 4) a *Gaussia* luciferase expression plasmid for normalization (pCMV-GLuc; 5 ng per well). After 24 h, cells were washed with 2x PBS solution and lysed with 75 µL of passive lysis buffer (Dual Luciferase Reporter Kit Assay; Promega) and then subjected to luciferase assay. *Firefly* luciferase activity (UAS-luc reporter constructs) was quantified using the Dual Luciferase Reporter Assay System (Promega) and normalized against cotransfected *Gaussia* luciferase expressing plasmid. Datasets are the results of at least four independent cell transfections and luciferase measurements. Statistical significance was assessed by means of the Wilcoxon-Mann-Whitney test (n = 4 or 8). Raw data are in Dataset S5.

**Data Availability.** The reference genome used in this study is deposited in European Nucleotide Archive (ENA) (PRJEB31270) and the sequence of the SDR together with its annotation as well as other sequences and codes used in this study are available on Github: https://github.com/nimarafati/Atlantic_herring_SDR. Variant calls used in the $χ^2$-analysis were from Ref. 19 deposited in Dryad (5r774). Also, Samplix and cDNA data generated for this study are deposited in ENA (PRJEB38031)

**ACKNOWLEDGMENTS.** We thank Carl-Henrik Heldin and Aristidis Moustakas for valuable comments and Chunheng Mo for technical assistance. The study was supported by the Knut and Alice Wallenberg Foundation (to L.A.), Vetenskapsrådet (to L.A.), Research Council of Norway project 254777, GENSINC (to A.F., F.B. and L.A.), and AquaCRISPR (ANR-16-COFA-0004-01) (to A.H.). The National Genomics Infrastructure (NGI)/Uppsala Genome Center and UPPMAX (Uppsala Multidisciplinary Center for Advanced Computational Science) provided service in massive parallel sequencing and computational infrastructure. Work performed at NGI/Uppsala Genome Center has been funded by RFI/VR (Council for Research Infrastructure/Vetenskapsrådet) and Science for Life Laboratory, Sweden.

1. D. Bachtrag et al.; Tree of Sex Consortium, Sex determination: Why so many ways of doing it? *PLoS Biol.* **12**, e1001899 (2014).

2. D. Charlesworth, J. E. Mank, The birds and the bees and the flowers and the trees: Lessons from genetic mapping of sex determination in plants and animals. *Genetics* **186**, 9–31 (2010).

3. G. S. van Doorn, Evolutionary transitions between sex-determining mechanisms: A review of theory. *Sex Dev.* **8**, 7–19 (2014).

4. A. Herpin, M. Schartl, Plasticity of gene-regulatory networks controlling sex determination: Of masters, slaves, usual suspects, newcomers, and usurpators. *EMBO Rep.* **16**, 1260–1274 (2015).

5. M. Schartl, D. Galiana-Arnoux, C. Schultheis, A. Böhne, J.-N. Volff, "A primer of sex determination in poeciliids" in *Ecology and Evolution of Poeciliid Fishes*, A. Pilastro, J. Evans, I. Schlupp, Eds. (The University of Chicago Press, 2010), pp. 264–275.

6. M. Matsuda, M. Sakaizumi, Evolution of the sex-determining gene in the teleostean genus Oryzias. *Gen. Comp. Endocrinol.* **239**, 80–88 (2016).

7. W. J. Gammerdinger, T. D. Kocher, Unusual diversity of sex chromosomes in African cichlid fishes. *Genes (Basel)* **9**, 480 (2018).

8. M. W. Pennell, J. E. Mank, C. L. Peichel, Transitions in sex determination and sex chromosomes across vertebrate species. *Mol. Ecol.* **27**, 3950–3963 (2018).

9. R. H. Devlin, Y. Nagahama, Sex determination and sex differentiation in fish: An overview of genetic, physiological, and environmental influences. *Aquaculture* **208**, 191–364 (2002).

10. M. Schartl, A. Herpin, "Sex determination in vertebrates" in *Encyclopedia of Reproduction*, M. K. Skinner, Ed., (Elsevier, 2018), pp. 159–167.

11. R. A. Fisher, The evolution of dominance. *Biol. Rev. Camb. Philos. Soc.* **6**, 345–368 (1931).

12. J. J. Bull, *Evolution of Sex Determining Mechanisms*, (Benjamin/Cummings, Menlo Park, CA, 1983).

13. W. R. Rice, The accumulation of sexually antagonistic genes as a selective agent promoting the evolution of reduced recombination between primitive sex-chromosomes. *Evolution* **41**, 911–914 (1987).

14. D. Charlesworth, B. Charlesworth, G. Marais, Steps in the evolution of heteromorphic sex chromosomes. *Heredity* **95**, 118–128 (2005).

15. M. Kirkpatrick, The evolution of genome structure by natural and sexual selection. *J. Hered.* **108**, 3–11 (2017).
16. H. Ida, N. Oka, K. Hayashigaki, Karyotypes and cellular DNA contents of three species of the subfamily Clupeinae. *J. Ichthyol.* **38**, 289-294 (1991).

17. M. E. Pettersson et al., A chromosome-level assembly of the Atlantic herring genome-detection of a supergene and other signals of selection. *Genome Res.* **29**, 1919-1928 (2019).

18. S. Lamichhaney et al., Parallel adaptive evolution of geographically distant herring populations on both sides of the North Atlantic Ocean. *Proc. Natl. Acad. Sci. U.S.A.* **114**, E3452-E3461 (2017).

19. A. Martinez Barrio et al., The genetic basis for ecological adaptation of the Atlantic herring revealed by genome sequencing. *eLife* 5, e12081 (2016).

20. E. B. Madsen, T. Kvist, I. Höijer, A. Ameur, M. J. Mikkelsen, Xdrop: Targeted sequencing of long DNA molecules from low input samples using droplet sorting bio-Rxiv. 10.1101/409086 (6 September 2018).

21. D. Ren et al., A sperm ion channel required for sperm motility and male fertility. *Nature* **413**, 603-609 (2001).

22. R. Seifert et al., The CatSper channel controls chemosensation in sea urchin sperm. *EMBO J.* **34**, 379-392 (2015).

23. C. Brenker et al., Action of steroids and plant triterpenoids on CatSper Ca²⁺ channels in human sperm. *Proc. Natl. Acad. Sci. U.S.A.* **115**, E344-E346 (2018).

24. J. F. Smith et al., Disruption of the principal, progesterone-activated sperm Ca²⁺ channel in a CatSper2-deficient infertile patient. *Proc. Natl. Acad. Sci. U.S.A.* **110**, 6823-6828 (2013).

25. M. S. Hildebrand et al., Genetic male infertility and mutation of CATSPER ion channels. *Eur. J. Hum. Genet.* **18**, 1178-1184 (2010).

26. Q. Pan et al., Identification of the master sex determining gene in Northern pike (*Esox lucius*) reveals restricted sex chromosome differentiation. *PLoS Genet.* **15**, e1008013 (2019).

27. A. Herpin, C. Lelong, P. Favrel, Transforming growth factor-β-related proteins: An ancestral and widespread superfamily of cytokines in metazoans. *Dev. Comp. Immunol.* **28**, 461-485 (2004).

28. S. Mawaribuchi, S. Yoshimoto, S. Ohashi, N. Takamatsu, M. Ito, Molecular evolution of vertebrate sex-determining genes. *Chromosome Res.* **20**, 139-151 (2012).

29. M. C. Adolfi, R. T. Nakajima, R. H. Nóbrega, M. Schartl, Intersex, hermaphroditism, and gonadal plasticity in vertebrates: Evolution of the Müllerian duct and Amh/Amhr2 signaling. *Annu. Rev. Anim. Biosci.* **7**, 149-172 (2019).

30. J. A. Visser, AMH signaling: From receptor to target gene. *Mol. Cell. Endocrinol.* **211**, 65-73 (2003).

31. L. Sèdes et al., Anti-Müllerian hormone recruits BMPR-IA in immature granulosa cells. *PLoS One* **8**, e81551 (2013).

32. X. Cai, D. E. Clapham, Evolutionary genomics reveals lineage-specific gene loss and rapid evolution of a sperm-specific ion channel complex: CatSpers and CatSperbeta. *PLoS One* **3**, e3569 (2008).

33. J.-J. Chung et al., CatSperζ regulates the structural continuity of sperm Ca²⁺ signaling domains and is required for normal fertility. *eLife* **6**, e23082 (2017).

34. D. Wachten, J. F. Jikeli, U. B. Kaupp, Sperm sensory signaling. *Cold Spring Harb. Perspect. Biol.* **9**, a028225 (2017).

35. U. B. Kaupp, T. Strünker, Signaling in sperm: More different than similar. *Trends Cell Biol.* **27**, 101-109 (2017).

36. P. V. Lishko et al., The control of male fertility by spermatozoan ion channels. *Annu. Rev. Physiol.* **74**, 453-475 (2012).

37. W. Bönigk et al., An atypical CNG channel activated by a single cGMP molecule controls sperm chemotaxis. *Sci. Signal.* **2**, ra68 (2009).

38. T. Strünker et al., A K⁺-selective cGMP-gated ion channel controls chemosensation of sperm. *Nat. Cell Biol.* **8**, 1149-1154 (2006).

39. C. Trötschel et al., Absolute proteomic quantification reveals design principles of sperm flagellar chemosensation. *EMBO J.* **39**, e102723 (2020).

40. R. Gauss, R. Seifert, U. B. Kaupp, Molecular identification of a hyperpolarization-activated channel in sea urchin sperm. *Nature* **393**, 583-587 (1998).

41. F. Windler et al., The solute carrier SLC9C1 is a Na⁺/H⁺-exchanger gated by an S4-type voltage-sensor and cyclic-nucleotide binding. *Nat. Commun.* **9**, 2809 (2018).

42. L. Wobig et al., A family of hyperpolarization-activated channels selective for protons. *Proc. Natl. Acad. Sci. U.S.A.* **117**, 13783-13791 (2020).

43. S. Fechner et al., A K(+)-selective CNG channel orchestrates Ca(2+)-signalling in zebrafish sperm. *eLife* 4, e07624 (2015).

44. A. Herpin et al., Transcriptional rewiring of the sex determining dmrt1 gene duplicate by transposable elements. *PLoS Genet.* **6**, e1000844 (2010).

45. C. L. Peichel et al., Assembly of the threespine stickleback Y chromosome reveals convergent signatures of sex chromosome evolution. *Genome Biol.* **21**, 177 (2020).

46. A. McKenna et al., The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. *Genome Res.* **20**, 1297-1303 (2010).

47. B. S. Weir, C. C. Cockerham, Estimating F-statistics for the analysis of population structure. *Evolution* **38**, 1358-1370 (1984).

48. P. Danecek et al.; 1000 Genomes Project Analysis Group, The variant call format and VCFtools. *Bioinformatics* **27**, 2156-2158 (2011).

49. B. L. Cantarel et al., MAKER: An easy-to-use annotation pipeline designed for emerging model organism genomes. *Genome Res.* **18**, 188-196 (2008).

50. M. Magrane; UniProt Consortium, UniProt knowledgebase: A hub of integrated protein data. *Database* (Oxford) **2011**, bar009 (2011).

51. M. Stanke, M. Diekhans, R. Baertsch, D. Haussler, Using native and syntenically mapped cDNA alignments to improve de novo gene finding. *Bioinformatics* **24**, 637-644 (2008).

52. C. Burge, S. Karlin, Prediction of complete gene structures in human genomic DNA. *J. Mol. Biol.* **268**, 78-94 (1997).

53. F. Sievers et al., Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. *Mol. Syst. Biol.* **7**, 539 (2011).

54. Team RC, *R: A Language and Environment for Statistical Computing*, (R Foundation for Statistical Computing, Vienna, Austria, 2019).

55. E. Paradis, K. Schliep, ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. *Bioinformatics* **35**, 526-528 (2019).

56. M. Suyama, D. Torrents, P. Bork, PAL2NAL: Robust conversion of protein sequence alignments into the corresponding codon alignments. *Nucleic Acids Res.* **34**, W609-W612 (2006).

57. D. L. Rabosky et al., An inverse latitudinal gradient in speciation rate for marine fishes. *Nature* **559**, 392-395 (2018).

58. Z. Yang, PAML 4: Phylogenetic analysis by maximum likelihood. *Mol. Biol. Evol.* **24**, 1586-1591 (2007).

59. M. Källberg et al., Template-based protein structure modeling using the RaptorX web server. *Nat. Protoc.* **7**, 1511-1522 (2012).

60. A. Chaikuad et al., Crystal structure of the cytoplasmic domain of the bone morphogenetic protein receptor type-1B (BMPR1B) in complex with FKBP12 and LDN-193189. Protein Data Bank. https://www.rcsb.org/structure/3MDY (accessed 3 February 2020).

61. C. A. Vines et al., Motility initiation in herring sperm is regulated by reverse sodium-calcium exchange. *Proc. Natl. Acad. Sci. U.S.A.* **99**, 2026-2031 (2002).

62. D. Raju et al., Accumulation of glucosylceramide in the absence of the beta-glucosidase GBA2 alters cytoskeletal dynamics. *PLoS Genet.* **11**, e1005063 (2015).

63. E. Birney, M. Clamp, R. Durbin, GeneWise and Genomewise. *Genome Res.* **14**, 988-995 (2004).

64. K. Katoh, K. Misawa, K. Kuma, T. Miyata, MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. *Nucleic Acids Res.* **30**, 3059-3066 (2002).

65. K. Katoh, D. M. Standley, MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. *Mol. Biol. Evol.* **30**, 772-780 (2013).

66. S. Guindon et al., New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. *Syst. Biol.* **59**, 307-321 (2010).

67. S. J. Grzegorski, E. F. Chiari, A. Robbins, P. E. Kish, A. Kahana, Natural variability of Kozak sequences correlates with function in a zebrafish model. *PLoS One* **9**, e108475 (2014).

68. Y. Hong et al., Establishment of a normal medakafish spermatogonial cell line capable of sperm production *in vitro*. *Proc. Natl. Acad. Sci. U.S.A.* **101**, 8011-8016 (2004).

69. F. Liu et al., A human Mad protein acting as a BMP-regulated transcriptional activator. *Nature* **381**, 620-623 (1996).

70. R. Betancur-R et al., *The Tree of Life and a New Classification of Bony Fishes*, (PLOS Currents Tree of Life, ed. 1, 2013).