The rise and fall of the ancient northern pike master sex-determining gene

To verify that the two copies of amh were not collapsed into a single sequence during assembly, we computed the total number of apparent heterozygous sites on the amh region in population genomic data with the expectation that the presence of two divergent gene copies should result in high apparent heterozygosity when remapped on the single copy assembled in the genome. In total, 106 variants were observed in the pool of E. lucius males on the amha region located between 12,906,561 bp and 12,909,640 bp on LG08 of E. lucius (GCA_004634155.1: CM015581.1), resulting from the mapping of reads originating from both amha and amhby, while only 12 variants (true allelic variations) were observed in the same region when mapping reads from the female pool originating only from amha. With male Pool-Seq reads from D. pectoralis, we observed only four variant sites on the ~3 kb amh region located between 23,447 bp and 25,910 bp on the flattened_line_2941 contig and zero variant sites from the female Pool-Seq reads. Compared to the ‘control’ of E. lucius Pool-Seq reads where one amh gene from female pool result in 12 variant sites and two amh gene from the male pool result in 106 variant sites, the low number of variant sites in both male and female pool of D. pectoralis support that only one amh gene is present in the genome of D. pectoralis regardless of the phenotypic sex of the individuals. Although sex-specific Pool-Seq reads were not available for U. pygmaea, we performed the same analysis with reads from the single male individual used to assemble the genome. No variant was observed in the ~3 kb region containing amh located between 760,763 bp and 757,962 bp on contig 633485 of our draft genome of U. pygmaea, supporting that only one amh gene is also present in this species.

To help estimate the size of the sex locus from sex-specific RAD markers, we determined the number of potential SbfI cleavage sites based on our draft genome assemblies for each species. For E. lucius (Canadian population), E. masquinongy (Iowa population), N. hubbsi, and D. pectoralis, we predicted the number of RAD-Seq cleavage sites present in each genome by counting the number of unambiguous matches for sequence of SbfI (CCTGCAGG), the restriction enzyme used in RAD-Seq library preparation (Herrera et al., 2015). On average, we found 31.8 RAD markers per Mb in E. lucius, 38.4 in E. masquinongy, 41.6 in N. hubbsi, 30.8 in D. pectoralis, and 26.2 in U. pygmaea. Because we do not see large species differences (26–40 RAD markers/Mb) this suggests that, apart from potential local variations of the RAD markers density in sex loci, our RAD-Seq comparative analysis could to some extent be used to compare sex locus size within species on the basis of this number of sex-specific RAD markers. We are aware of the limitation that using the number of sex-specific marker usually lead to an overestimation of the size of the sex locus. For all of our species with the exception of U. pygmaea, we identified very few sex-specific markers, indicating very small sex locus. This simple calculation is only intended to helps provide a rough approximation of the size of the sex locus.

*Samples from (Ouellet-Cauchon et al., 2014). **Sex was recorded in this E. masquinongy population based on the urogenital pores morphology. ***Sex was recorded in N. hubbsi based on the specific coloration of males during the breeding season. NR = phenotypes not recorded. NA = not applicable (sex phenotypes not recorded). Museum collection numbers are as follows: MNHN 2014–2719, MNHN 2014–2720, MNHN 2014–2721, MNHN 2014–2722, and MNHN 2014–2723 for E. cisalpinus and MNHN 2013–1246, MNHN 2013–1245, and MNHN 2013–838 for E. aquitanicus.

The number of markers retained correspond to the number of markers present with depth higher than min. depth in at least one individual.