The rise and fall of the ancient northern pike master sex-determining gene
Abstract
The understanding of the evolution of variable sex determination mechanisms across taxa requires comparative studies among closely related species. Following the fate of a known master sex-determining gene, we traced the evolution of sex determination in an entire teleost order (Esociformes). We discovered that the northern pike (Esox lucius) master sex-determining gene originated from a 65 to 90 million-year-old gene duplication event and that it remained sex linked on undifferentiated sex chromosomes for at least 56 million years in multiple species. We identified several independent species- or population-specific sex determination transitions, including a recent loss of a Y chromosome. These findings highlight the diversity of evolutionary fates of master sex-determining genes and the importance of population demographic history in sex determination studies. We hypothesize that occasional sex reversals and genetic bottlenecks provide a non-adaptive explanation for sex determination transitions.
Introduction
Genetic sex determination (GSD) evolved independently and repeatedly in diverse taxa, including animals, plants, and fungi (Tree of Sex Consortium et al., 2014; Beukeboom and Perrin, 2014), but the stability of such systems varies drastically among groups. In mammals and birds, conserved male or female heterogametic sex determination (SD) systems have been maintained over a long evolutionary time with conserved master sex-determining (MSD) genes (Marshall Graves, 2008). In contrast, teleost fishes display both genetic and environmental sex determination (ESD), and a remarkable evolutionary lability driven by rapid turnovers of sex chromosomes and MSD genes (Kikuchi and Hamaguchi, 2013; Pan et al., 2017). These characteristics make teleosts an attractive group in which to study the evolution of SD systems.
In the past two decades, the identity of a variety of MSD genes has been revealed in teleosts thanks to advances in sequencing technologies. These findings generated new hypotheses on how the birth of MSD genes through allelic diversification or duplication with or without translocation may drive sex chromosome turnover in vertebrates (Kikuchi and Hamaguchi, 2013; Pan et al., 2017). These teleost MSD genes also provided empirical support for the ‘limited option’ idea that states that certain genes known to be implicated in sex differentiation pathways are more likely to be recruited as new MSD genes (Marshall Graves and Peichel, 2010). The majority of these recently discovered MSD genes, however, were phylogenetically scattered, making it challenging to infer evolutionary patterns and conserved themes of sex chromosomes and/or MSD gene turnovers. Although comparative studies have been accomplished in medakas, poeciliids, tilapiine cichlids, salmonids, and sticklebacks (Kikuchi and Hamaguchi, 2013), transitions of SD systems in relation to the fate of known MSD genes within closely related species have only been explored in medakas (Myosho et al., 2015) and salmonids (Guiguen et al., 2018).
Esociformes is a small order of teleost fishes (Figure 1) that diverged from their sister group Salmoniformes about 110 million years ago (Mya) and diversified from a common ancestor around 90 Mya (Campbell et al., 2013; Campbell and Lopéz, 2014). With two families (Esocidae, including Esox, Novumbra, and Dallia, and Umbridae, including Umbra) and 13 well-recognized species (Warren et al., 2020), Esociformes is an ecologically important group of freshwater species from the northern hemisphere (Campbell et al., 2013; Campbell and Lopéz, 2014). We demonstrated that a male-specific duplication of the anti-Müllerian hormone gene (amhby) is the MSD gene in Esociformes and that this gene is located in a small male-specific insertion on the Y chromosome of northern pike (Esox lucius) (Pan et al., 2019).
Here, we took advantage of the small number of Esociformes species and their relatively long evolutionary history to explore the evolution of SD in relation to the fate of this amhby MSD gene. Generating novel draft genome assemblies and population genomic data in multiple species of Esociformes (Figure 1), we traced the evolutionary trajectory of amhby from its origin in a gene duplication event 65–90 Mya to its demise during species- or population-specific SD transitions. Our results highlight the diverse evolutionary fates experienced by an MSD gene. Among SD transitions, we uncovered a recent loss of the entire Y chromosome in some populations of northern pike. We hypothesize that drift, exacerbated by the bottleneck effect, may have facilitated the loss of the sex chromosome along with its MSD gene in these populations, a new mechanism for entire sex chromosome loss in vertebrates.
Results
The Esox lucius amhby gene originated from an ancient duplication event
Previously, we demonstrated that a male-specific duplication of the anti-Müllerian hormone gene (amhby) is the MSD gene in Esociformes and that this gene is located in a small, ~140 kb male-specific insertion on the Y chromosome of northern pike (Esox lucius) (Pan et al., 2019). To explore the evolution of the autosomal copy of amh (amha) and amhby in Esociformes, we collected phenotypically sexed samples of most species of this order (Figure 1 and Table 1). We searched for homologs of amh using homology cloning and whole-genome sequencing. We found two amh genes in all surveyed Esox and Novumbra species. In the more basally diverging species, Dallia pectoralis and Umbra pygmaea, we found only one amh gene in both species in tissue-specific transcriptome databases (Pasquier et al., 2016), while two amh transcripts were readily identified in E. lucius (Pan et al., 2019).
To clarify relationships among these amh homologs, we inferred phylogenetic trees from these sequences. These phylogenies provided a clear and consistent topology with two well-supported gene clusters among the amh sequences from Esocidae (Figure 2A, Appendix 1—figure 1). The first gene cluster contains the autosomal amha of E. lucius and an amh homolog from all other species of Esox, Dallia, and Novumbra. Therefore, we assigned these amh homologs as amha orthologs. The second gene cluster contains the Y-specific amhby gene of E. lucius, along with the other amh homolog from all species for which we identified two amh sequences. We assigned these amh homologs as amhby orthologs. The general topology of each cluster was in agreement with the Esociformes species taxonomy.
In U. pygmaea, the closest sister species to the Esocidae, we found a single amh homolog, which roots at the base of the amha/amhby clusters. This result indicates that the amh duplication happened after the divergence of Esocidae (including Esox, Dallia, and Novumbra) and Umbridae (including Umbra) lineages. In contrast to other Esocidae that have two amh orthologs, we found only one amh gene in D. pectoralis, which fell in the phylogeny within the amha ortholog cluster, suggesting that this species experienced a secondary loss of its amhby gene after the amh duplication at the root of the Esocidae lineage (Figure 2B and Appendix 1—figure 1). We confirmed the absence of an additional divergent amh gene in D. pectoralis by searching sex-specific Pool-Seq reads from 30 males and 30 females. In addition, only one amh homolog was found in an ongoing genome assembly project with long reads for a male D. pectoralis (Y. Guiguen, personal communication). Using information from an Esociformes time-calibrated phylogeny (Campbell et al., 2013), we estimated that the amhby duplicate arose between 65 and 90 Mya (Figure 2B).
Sex linkage of amhby in the Esocidae lineage
Given that amhby originated from an ancient duplication event and is the MSD gene in E. lucius, we investigated when amhby became the MSD gene and whether transitions of SD systems exist in this group. Because male sex linkage of a sequence is a strong indication that this sequence is located on a Y-chromosome-specific region, we first investigated the association of amhby with male phenotype. We found an association between male phenotype and the presence of amhby in most species of Esox and Novumbra (Table 1). This sex linkage was significant in European, Asian, and Alaskan populations of E. lucius, in two North American populations of E. masquinongy (Iowa and Wisconsin), and in all populations surveyed of E. reichertii, E. americanus, E. niger, and N. hubbsi. For two recently described species, E. cisalpinus and E. aquitanicus (Denys et al., 2014), we had insufficient samples with clear species and sex identification for a decisive result. We confirmed the presence of amhby in all males and its absence from all females in these two species, but the association was not significant due to low sample size (Table 1). Because amhby is associated with male phenotype in most species of Esox and Novumbra, it likely gained an MSD role shortly after its origin either at the root of Esocidae or before the split of Esox and Novumbra lineages (~56 Mya) (Campbell et al., 2013). Despite this global conservation of the linkage between amhby and male phenotype, however, we found some variations in amhby sex association across populations of E. lucius and E. masquinongy (Table 1). These population variations are further investigated below.
Whole-genome analyses of the evolution of sex determination systems in Esociformes
Because amhby was not completely associated with male phenotype in E. masquinongy and N. hubbsi (Table 1) and the gene was not found in D. pectoralis and U. pygmaea, we also used population genomic approaches to search for whole-genome sex-specific signatures in these species. Because of the close phylogenetic distance (45 Mya) between E. lucius and E. masquinongy, we used the E. lucius genome assembly to remap reads from E. masquinongy. We performed Pool-Seq analysis of E. masquinongy (Iowa, USA) to compare the size and location of its sex locus with that of E. lucius (Pan et al., 2019). Comparison of the sex-specific heterozygosity across the entire genome revealed that a single region of less than 50 kb, containing the highest density of male-specific single-nucleotide polymorphisms (SNPs) in E. masquinongy, is located in a region homologous to the proximal end of LG24, where the sex locus of E. lucius is located. Besides LG24, no other linkage group showed enrichment for sex-biased chromosome differentiation (Appendix 1—figure 2), suggesting that the location and the small size of the sex locus are likely conserved between E. lucius and E. masquinongy.
No high-quality genome assembly was available from a closely related species to N. hubbsi. Therefore, we performed de novo RAD-Seq analysis on 21 males and 19 females for N. hubbsi. We found two markers showing a significant association with male sex and no female-biased marker (Figure 3A). This result supports the male heterogametic SD system (XX/XY) that was suggested by the amhby sex linkage and also indicates that the N. hubbsi sex locus is likely small (41.6 restriction enzyme cutting sites/Mb across the genome, on average, Supplementary file 2).
-
Figure 3—source code 1
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-code1-v2.zip
-
Figure 3—source data 1
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-data1-v2.csv
-
Figure 3—source data 2
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-data2-v2.csv
-
Figure 3—source data 3
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-data3-v2.csv
-
Figure 3—source data 4
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-data4-v2.csv
-
Figure 3—source data 5
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig3-data5-v2.csv
In D. pectoralis and U. pygmaea, two species in which amhby is absent, we carried out RAD-Seq analyses comparing phenotypic males and females. In D. pectoralis, we found only three female-biased RAD markers, suggesting that this species also has a small sex locus region (30.8 restriction enzyme cutting sites/Mb on average, Supplementary file 2) under a female heterogametic SD system (ZZ/ZW) (Figure 3B). This ZZ/ZW SD system was further supported by an independent Pool-Seq analysis revealing 45 times more female-specific k-mers (N = 1,081,792) than male-specific k-mers (N = 23,816). This excess of female-specific k-mers indicates that females carry genomic regions that are absent from males and thus that females are the heterogametic sex. In contrast, 140 male-biased and no female-biased RAD markers were identified in U. pygmaea, supporting that this species has a large sex locus (26.2 restriction enzyme cutting sites/Mb on average, Supplementary file 2) under a male heterogametic SD system (XX/XY) (Figure 3C).
Some populations of Esox lucius lost their Y chromosome and ancestral master sex-determining gene
Although amhby is the MSD gene in European populations of northern pike (Pan et al., 2019), this gene was absent from the genome assembly (GCA_000721915.3) of a male specimen from a Canadian population (GCA_000721915.3). To explore this discrepancy, we surveyed the sex linkage of amhby in geographically isolated populations of E. lucius. We found significant male linkage in all investigated European and Asian populations and in one Alaskan population (Table 1). In contrast, amhby was absent from both males and females in all other North American populations.
To investigate whether the loss of amhby in most North American populations coincides with large genomic changes, we compared genome-wide sex divergence patterns between a European population carrying the amhby gene (Ille-et-Vilaine, France) and a North American population that lacks amhby (Quebec, Canada) using a Pool-Seq approach. We aligned Pool-Seq reads to an improved European E. lucius genome assembly (NCBI accession number SDAW00000000; assembly metrics presented in Supplementary file 3) in which all previously identified Y-specific contigs (Pan et al., 2019) were scaffolded into a single contiguous locus on the Y chromosome (LG24). In the European population, Pool-Seq analysis confirmed, with better resolution, previous results (Pan et al., 2019), showing the presence of a small Y chromosome region (~140 kb) characterized by many male-specific sequences at the proximal end of LG24 (Figure 4 A.1 and 4B.1). In comparison, virtually no reads from either male or female pools from the Quebec population mapped to this 140 kb Y-specific region (Figure 4 A.2 and 4B.2). Furthermore, we observed no differentiation between males and females along the remainder of the genome with either Pool-Seq or reference-free RAD-Seq (Appendix 1—figure 3). Together, these results suggest that the Quebec population and likely other mainland US populations, where the MSD gene amhby was not found, lack not only amhby, but also the surrounding Y-specific region identified in European populations. Moreover, the new sex locus of these North American populations, if it exists, is too small to be detected by the RAD-Seq and Pool-Seq approaches.
-
Figure 4—source code 1
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-code1-v2.zip
-
Figure 4—source data 1
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-data1-v2.csv
-
Figure 4—source data 2
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-data2-v2.csv
-
Figure 4—source data 3
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-data3-v2.csv
-
Figure 4—source data 4
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-data4-v2.csv
-
Figure 4—source data 5
- https://cdn.elifesciences.org/articles/62858/elife-62858-fig4-data5-v2.csv
Evolving sex determination systems: the case of the muskellunge, Esox masquinongy
As in E. lucius, we found variation in amhby sex linkage among different populations of E. masquinongy (Table 1). In addition, we found that two males and one female of E. masquinongy were heterozygous for amhby in the population from Quebec (Appendix 1—figure 4), a finding that conflicts with the expected hemizygous status of a Y-specific MSD gene. We thus used RAD-Seq to compare genome-wide patterns of sex differentiation in a population from Iowa (USA) where amhby was significantly associated with male phenotype and in a population from Quebec (Canada) where it was not. In the Iowa population, five RAD markers were significantly associated with male phenotype, while no marker was associated with female phenotype (Figure 3D). This result indicates a male heterogametic SD system (XX/XY) with a low differentiation between the X and Y chromosomes, as observed in northern pike (Pan et al., 2019) and N. hubbsi (this study). In contrast, we did not find a sex-specific marker in the Quebec population (Figure 3E). This result was further supported by the analysis of a publicly available dataset from another Quebec population (PRJNA512459, Appendix 1—figure 5), suggesting that the sex locus of E. masquinongy from these Quebec populations is too small to be detected with RAD-Seq (38.4 restriction enzyme cutting sites/Mb on average, Supplementary file 2) or that this population displays a multifactorial GSD or ESD system.
Evolution of the structure of the amhby gene in the Esocidae
The typical amh gene in teleost fishes comprises seven exons encoding a protein that contains 500–571 amino acids. The conserved C-terminal TGF-β domain is crucial for canonical Amh function (di Clemente et al., 2010). The predicted structures of most of the amha and amhby genes in Esociformes are consistent with this canonical structure, and most amha and amhby genes do not show signature of relaxation from purifying selection (Supplementary file 4 and Supplementary file 5). However, in both E. niger and N. hubbsi, the predicted Amhby protein is truncated in its N- or C-terminal part. In E. niger, this truncated amhby gene is flanked by repeated elements and encompasses only three of the seven conserved Amh exons, that is exons 5, 6, and 7 with a truncated TGF-β domain (Appendix 1—figure 6). In N. hubbsi, the truncated amhby gene contains the seven conserved Amh exons but with an N-terminal truncation of exon 1 encoding only eight amino acids with no homology to the conserved 50 amino acid sequence of the first Amh exon of other Esocidae (Appendix 1—figure 6). Together, these results show that in some species where amhby is sex linked, the protein sequence is strongly modified in the N- and C-terminal-conserved domains that are needed for proper protein secretion and correct conformation (di Clemente et al., 2010), leaving questions as to whether these proteins are functional.
Discussion
The northern pike MSD gene amhby substantially diverges from its autosomal copy amha, suggesting an ancient origin (Pan et al., 2019), unlike other cases of amh duplication in fishes, that appears to be comparatively young (Hattori et al., 2013; Li et al., 2015). Here, we show that the amhby duplicate emerged before the split of Esocidae and Umbridae between 65 and 90 Mya (Campbell et al., 2013) and was subsequently secondarily lost in the Dallia lineage. We also demonstrate that amhby presence is significantly associated with male phenotype in most species of Esox and in Novumbra. Although we cannot rule out the possibility that amhby was recruited independently as the MSD gene in Esox and Novumbra, our results suggest the more parsimonious hypothesis that amhby likely acquired an MSD function before the diversification of Esocidae at least 56 Mya.
Studies on SD systems of cichlid fishes and true frogs suggested that frequent turnovers can keep sex chromosomes undifferentiated (Gammerdinger et al., 2018; Jeffries et al., 2018). It was proposed that in these organisms, the turnover of sex chromosomes was facilitated by autosomes that host genes involved in sex differentiation pathways that could be readily recruited as new MSD genes, facilitating turnovers of sex chromosomes (Vicoso, 2019). Alternative evolutionary pathways could also involve conserved MSD genes that translocate onto different autosomes, as demonstrated for the salmonid sdY gene (Guiguen et al., 2018). Interestingly, this scenario does not seem to be the case with Esocidae. No highly differentiated sex chromosomes were observed in the surveyed population of E. masquinongy and N. hubbsi (this study), or in E. lucius, despite Esocidae having retained the same ancestral XX/XY system (likely controlled by amhby) for at least 56 million years. In addition, whole-genome analyses in E. masquinongy suggest that its male-specific locus on the Y chromosome is conserved in terms of size and location with that of E. lucius despite 45 million years of divergence. Such conservation is unusual in teleosts, where frequent turnover of SD systems is considered the norm (Kikuchi and Hamaguchi, 2013; Pan et al., 2017). In line with our present results, substantially undifferentiated sex chromosomes have also been maintained over relatively long evolutionary periods in some Takifugu fish species (Kamiya et al., 2012), supporting the idea that theoretical models of sex chromosome evolution cannot be generalized across taxa and that additional, as-yet unknown evolutionary forces can prevent sex chromosome decay.
Even though this MSD gene has been maintained for a long evolutionary period, our results revealed that it has been subjected to several independent SD transitions, potentially driven by different mechanisms in the Esocidae lineage. In D. pectoralis, we observed a transition from the ancestral XX/XY male heterogametic system shared by Esocidae and U. pygmaea to a ZZ/ZW female heterogametic system. This transition also coincided with the loss of amhby. This shift from male to female heterogamety could have been driven by the rise of a new dominant female determinant, leading to the obsolescence and the subsequent loss of the ancestral male MSD gene, similar to the transition of SD documented among the Oryzias and also in the housefly (Kozielska et al., 2008, Myosho et al., 2012).
The other SD transitions that we characterized, those in E. lucius and E. masquinongy, are probably more recent because these transitions are found within populations of the same species. In E. masquinongy, we observed population-specific variation in SD systems, supporting the conservation of an XX/XY system with amhby as the MSD gene in two US populations belonging to the northern lineage, but not in a Quebec population belonging to the Great Lakes lineage (Turnquist et al., 2017). This result suggests that although amhby is still present in the genome of the Great Lakes lineage, it is no longer the MSD gene and that a new population-specific SD mechanism emerged after the split of the lineages. Previously, female heterogamety was inferred for the Great Lakes lineage based on the presence of males in the offspring of gynogenetic females (Dabrowski et al., 2000). With RAD-Seq, we did not find evidence for female-specific genomic regions in the Great Lakes lineage of E. masquinongy. Higher resolution approaches are needed to identify differentiation between male and female genomes to reveal the new SD locus and its interaction with amhby in the Great Lakes lineage.
In E. lucius, we found that amhby is present and completely sex linked in European, Asian, and Alaskan regions, but is not present elsewhere in North America (Figure 5). This finding suggests that amhby functions as an MSD gene in Eurasian and Alaskan populations, but that populations from elsewhere in North America lost amhby together with the entire sex locus. The current geographic distribution of E. lucius results from a post-glacial expansion from three glacial refugia ~0.18 to 0.26 Mya. The North American populations lacking amhby belong to a monophyletic lineage with a circumpolar distribution originating from the same refugia as the other population we surveyed that carry amhby (Skog et al., 2014). Fossil records from Alaska support the idea of an ‘out of Alaska’ North American expansion of E. lucius within the last 100,000 years (Wooller et al., 2015), suggesting that the Y-specific sequences containing amhby were lost during this dispersal period (Figure 5, Appendix 1—figure 7).
The loss of the entire sex locus in such a short evolutionary time is unlikely to have resulted from the pseudogenization of the amhby gene, as is the case for the ancestral MSD gene of the Luzon medaka (Myosho et al., 2012). Rather, it is probably the result of colonization by a small pool of females and sex-reversed XX males. Such a founder-effect hypothesis is supported by the fact that North American populations of E. lucius display a much lower genetic diversity compared to other populations from the same lineage (Skog et al., 2014). This hypothesis is also supported by the fact that environmental influence on GSD systems is a well-documented phenomenon in fishes (Baroiller et al., 2009; Chen et al., 2014; Sato et al., 2005) and that XX sex reversal induced by environmental factors has been observed occasionally in captive E. lucius (Pan et al., 2019). Our whole-genome analyses failed to identify any sex-associated signals in these North American populations, meaning that if such a new sex locus exists, it would lack detectable signatures of sequence differentiation, similar to the one-SNP sex locus of some Takifugu species (Kamiya et al., 2012). Whether ESD, which is common among poikilothermic animals including teleosts (Navara, 2018), could be the SD system in these North American populations or whether these populations already have a new GSD system in place is still unresolved. Notably, this Y chromosome loss in natural populations of E. lucius mirrors the loss of the W chromosome in lab strains of zebrafish after just a few decades of selective breeding (Wilson et al., 2014).
A few evolutionary models have been formulated to capture the dynamics of sex chromosome turnover (Vicoso, 2019). The replacement of the ancestral SD locus by a new one could be driven by drift alone, by positive selection (Bull and Charnov, 1977), by sexual antagonistic selection (van Doorn and Kirkpatrick, 2007), or by the accumulation of deleterious mutations at the sex locus (Blaser et al., 2013). In all of these models, the ancestral sex chromosomes would revert to autosomes only when the new SD locus is fixed in the population. However, in North American populations of E. lucius, the ancestral sex chromosome was likely lost due to drift in bottlenecked populations. Interestingly, this process would not require the simultaneous emergence of a new GSD system, given the flexibility of SD mechanisms in teleosts, where environmental cues generate phenotypic sexes in the absence of a GSD system. ESD, which is easily invaded by a new genetic system (Muralidhar and Veller, 2018; van Doorn, 2014), may serve as a transitional state between sex chromosome turnovers.
We found two cases of altered amhby copies. In N. hubbsi and E. niger, amhby is strongly sex linked, which suggests that it is located in a Y-specific region of the genome. In both species, however, the amhby gene lacks parts of the conserved C- and/or N-terminal Amh regions, which are needed for proper protein secretion and interaction between Amh and its receptor (di Clemente et al., 2010). Although we cannot exclude that these are non-functional amhby copies, this hypothesis seems unlikely because the tight sex linkage of amhby in these species would require the independent emergence of a new MSD gene in close proximity to amhby. It is more likely that these copies still function as MSD genes, but with altered protein sequences. A similar example of a truncated amh duplicate acting as an MSD gene is found in the teleost cobalt silverside Hypoatherina tsurugae (Bej et al., 2017). Additional cases of truncated duplicated genes functioning as MSD genes have been described in Salmonids, yellow perch (Perca flavescens), and Atlantic herring (Clupea harengus), suggesting that the preservation of ancestral domains is not necessary for a duplicated protein to assume a novel role (Feron et al., 2020b; Rafati et al., 2020; Yano et al., 2012). These cases are all consistent with domain gains and losses contributing to new protein functions (Moore and Bornberg-Bauer, 2012).
Collectively, our results depict the evolutionary trajectories of a conserved MSD gene in an entire order of vertebrates, highlighting both the potential stability of MSD genes as well as non-differentiated sex chromosomes in some lineages, and the dynamics of species- or population-specific SD evolution in teleost fishes. Our results highlight the importance of careful consideration of the population demographic history of SD systems, and of the potential buffering role of ESD during transitions between genetic SD systems in sex chromosome evolution.
Materials and methods
Sample collection
Request a detailed protocolInformation on the Esociformes species, collectors, sexing method, and experiments is provided in Supplementary file 6.
Genomic DNA extraction
Request a detailed protocolFin clips were preserved in 75–100% ethanol at 4°C until genomic DNA (gDNA) extraction. For genotyping, samples were lysed with 5% Chelex and 25 mg of Proteinase K at 55°C for 2 hr, followed by incubation at 99°C for 10 min. For Illumina sequencing, gDNA was obtained using NucleoSpin Kits for Tissue (Macherey-Nagel, Duren, Germany) following the producer’s protocol. DNA concentration was quantified using Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA) and a Qubit3 fluorometer (Invitrogen, Carlsbad, CA). For Pool-Seq, DNA from different samples was normalized to the same quantity before pooling for male and female libraries separately. High-molecular-weight gDNA for long-read sequencing was extracted as described by Pan et al., 2019.
Genome and population genomics sequencing
Request a detailed protocolDraft genomes of northern pike, muskellunge (E. masquinongy), chain pickerel (E. niger), Olympic mudminnow (N. hubbsi), and Alaska blackfish (D. pectoralis) were sequenced using a whole-genome shotgun strategy with 2 × 250 bp Illumina reads. Libraries were built using the Truseq nano kit (Illumina, FC-121–4001) following instructions from the manufacturer. gDNA (200 ng) was briefly sonicated using a Bioruptor (Diagenode). The gDNA was end-repaired and size-selected on beads to retain fragments of ~550 bp, and these fragments were a-tailed and ligated to Illumina’s adapter. The ligated gDNA was then subjected to eight PCR cycles. Libraries were checked on a fragment analyzer (AATI) and quantified by qPCR using a library quantification kit from KAPA. Libraries were sequenced on a HIseq2500 using the paired end 2 × 250 bp v2 rapid mode according to the manufacturer’s instructions. Image analysis was performed by HiSeq Control Software, and base calling was achieved using Illumina RTA software. The output of each run is presented in Supplementary file 7. To improve the European male genome of E. lucius, we generated an extra coverage of Oxford Nanopore long reads using a higher fragment size (50 kb) library made from gDNA extracted from a different male from the same European population as the sample used for the previous genome assembly. Library construction and genome sequencing were carried out as in Pan et al., 2019 and 12.7 Gbp of new data were generated from one PromethION flowcell.
Pool-Seq was performed on E. masquinongy, D. pectoralis, and the North American populations of E. lucius. Pooled libraries were constructed using a Truseq nano kit (Illumina, FC-121–4001) following the manufacturer’s instructions. Male and female DNA Pool-Seq libraries were prepared for each species using the Illumina TruSeq Nano DNA HT Library Prep Kit (Illumina, San Diego, CA) with the same protocol as for the draft genome sequencing. The libraries were then sequenced on a NovaSeq S4 lane (Illumina, San Diego, CA) using the paired-end 2 × 150 bp mode with Illumina NovaSeq Reagent Kits following the manufacturer’s instructions. The output of each run for each sex is presented in Supplementary file 7.
RAD-Seq was performed on E. masquinongy, N. hubbsi, D. pectoralis, U. pygmaea, and the North American populations of E. lucius. RAD libraries were constructed from gDNA extracted from fin clips for each species using a single Sbf1 restriction enzyme as in Amores et al., 2011. Each library was sequenced on one lane of an Illumina HiSeq 2500. The summary of the output of each dataset is presented in Supplementary file 8.
Analysis of population genomic data for sex-specific signals
Request a detailed protocolRaw RAD-seq reads were demultiplexed with the process_radtags.pl script from stacks version 1.44 (Catchen et al., 2013) with all parameters set to default. Demultiplexed reads for each species were analyzed with the RADSex computational workflow (Feron et al., 2020a) using the radsex software version 1.1.2 (10.5281/zenodo.3775206). After generating a table of markers depth with process, the distribution of markers between sexes was computed with distrib and markers significantly associated with phenotypic sex were extracted with signif using a minimum depth of 10 (−d 10) for both commands and with all other settings at the default. RAD-Seq figures were generated with sgtr (10.5281/zenodo.3773063).
For each Pool-Seq dataset, reads were aligned to the reference genome using BWA mem (version 0.7.17) (Li and Durbin, 2009). Alignment results were sorted by genomic coordinates using samtools sort (version 1.10) (Li et al., 2009), and PCR duplicates were removed using samtools rmdup. A file containing nucleotide counts for each genomic position was generated with the pileup command from PSASS (version 3.0.1b, 10.5281/zenodo.3702337). This file was used as input to compute FST between males and females, number of male- and female-specific SNPs, and male and female depth in a sliding window along the entire genome using the analyze command from PSASS with parameters: window-size 50,000, output-resolution 1000, freq-het 0.5, range-het 0.15, freq-hom 1, range-hom 0.05, min-depth 1, and group-snps. The analysis was performed with a snakemake workflow available at https://github.com/SexGenomicsToolkit/PSASS-workflow. Pool-Seq figures were generated with sgtr version 1.1.1 (10.5281/zenodo.3773063).
For k-mer analysis, 31-mers were identified and counted in the reads of the male and female pools using the ‘count’ command from Jellyfish (version 2.2.10) (Marçais and Kingsford, 2011) with the option ‘-C’ activated to count only canonical k-mers and retaining only k-mers with an occurrence >5 and <50,000,000. Tables of k-mer counts produced by Jellyfish were merged with the ‘merge’ command from Kpool (https://github.com/INRA-LPGP/kpool), and the resulting merged table was filtered using the ‘filter’ command to only retain k-mers present >25 times in one sex and <5 times in the other sex; such k-mers were considered sex biased.
Sequencing of amha and amhby genes
Request a detailed protocolTo survey the presence of amha (the canonical copy of amh) and amhby in Esociformes, we collected samples of 11 species (Table 1): all species in the genus Esox (E. americanus americanus, E. americanus vermiculatus, E. aquitanicus, E. cisalpinus, E. masquinongy, E. reichertii, and E. niger), N. hubbsi (the only species in its genus), D. pectoralis (the only well-recognized species in its genus), and one species from the genus Umbra (U. pygmaea). To search for amh homologs in the genomes of these 11 species, we either sequenced PCR amplicons with custom primers (Supplementary file 9) and/or sequenced and assembled the genome (Supplementary file 10) and searched in the assembly and in the raw reads for the presence of one or two amh genes.
For species closely related to E. lucius (E. aquitanicus, E. cisalpinus, and E. reichertii), amh homologs were amplified with primers designed on amha or amhby of E. lucius. Complete sequences of amh homologs were obtained from overlapping amplicons covering the entire genomic regions of both amha and amhby with primers anchored from upstream of the start codon (SeqAMH2Fw1 and SeqAMH1Fw1) and downstream of the stop codon (SeqAMH2Rev3 and SeqAMH1Rev4). All PCR amplicons were then Sanger sequenced from both directions and assembled to make consensus gene sequences.
For three other Esocidae species (E. niger, E. masquinongy, N. hubbsi), only partial amplifications of amh sequences were obtained with primers designed on amha or amhby of E. lucius. To acquire the complete amh sequences, we generated draft genome assemblies from amhby-carrying individuals.
For the two species from the more divergent genera (D. pectoralis and U. pygmaea), we were unable to amplify amhby with primers designed on regions that appear conserved on Esox species. Therefore, we also generated a genome assembly from a phenotypic male of each species.
For E. niger, E. masquinongy, N. hubbsi, and D. pectoralis, we used amha and amhby genomic sequences from E. lucius as queries and searched for amh homologs with blastn (https://blast.ncbi.nlm.nih.gov/Blast.cgi, version 2.10.0+) (Altschul et al., 1997) using the parameters ‘Max target sequences’=50 and ‘Max matches in a query range’=1.
For U. pygmaea, the most divergent species from E. lucius in this study, blasting with E. lucius amha and amhby sequences did not yield any result. We used the blastn strategy with the coding sequence of amh as well as tblastn of the protein sequence from Salmon salar, which returned only one contig.
No genome was available for the other two more distantly related species of E. lucius, E. americanus americanus, and E. americanus vermiculatus; therefore, we designed primers on regions conserved in multi-sequence alignments of amha and amhby in the other Esox species (Supplementary file 9).
To investigate whether truncations of Amhby in E. niger and N. hubbsi could be assembly artifacts, we searched for potential missing homologous sequences in raw genome reads of both species using tblastn (https://blast.ncbi.nlm.nih.gov/Blast.cgi, version 2.10.0+) (Altschul et al., 1997) using the exon 1 protein sequence from E. lucius as a query. Homologous sequences were determined with a maximum e-value threshold of 1E-10.
Validating amh gene number in D. pectoralis and U. pygmaea
Request a detailed protocolTo check whether the two copies of amh were not collapsed into a single sequence during assembly, we computed the total number of apparent heterozygous sites in 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. We compared these heterozygosity data in species where we only identified one amh gene, that is D. pectoralis and U. pygmaea, to results observed when mapping sex-specific pooled libraries of E. lucius (n = 30 males and 30 females, respectively) to a female reference genome containing only amha (GCA_004634155.1). Because there are two amh genes in male genomes of E. lucius and one amh gene in female genomes, this result from mapping E. lucius sex-specific Pool-Seq reads serves as a reference for expected number of variant sites. Reads from the male and female pools of E. lucius and from male pool of D. pectoralis were aligned separately to the reference genome (GCA_004634155.1) and our draft genome, respectively, using BWA mem version 0.7.17 (Li and Durbin, 2009) with default parameters. Each resulting BAM file was sorted and PCR duplicates were removed using Picard tools version 2.18.2 (http://broadinstitute.github.io/picard) with default parameters. Variants were then called for each BAM file using bcftools mpileup version 1.9 (Li et al., 2009) with default parameters on genomic regions containing amh, which locates between 12,906,561 and 12,909,640 bp on LG08 of E. lucius (GCA_004634155.1: CM015581.1), between 23,447 and 25,910 bp on the flattened_line_2941 contig of the draft genome of D. pectoralis, and between 760,763 and 757,962 bp on contig 633485 of our draft genome of U. pygmaea.
Gene phylogenetic analysis
Request a detailed protocolPhylogenetic reconstructions were performed on all amh homolog sequences obtained from Esociformes with the S. salar amh used as an outgroup. Full-length CDS were predicted with the FGENESH+ (Solovyev et al., 2006) suite based on the genomic sequence and Amh protein sequence from E. lucius. Sequence alignments were performed with MAFFT (version 7.450) (Katoh and Standley, 2013). Both maximum-likelihood and Bayesian methods were used for tree construction with IQ-TREE (version 1.6.7) (Minh et al., 2020) and Phylobayes (version 4.1) (Lartillot et al., 2007), respectively.
Additional analyses were carried with or without the truncated amh/Amh sequences of E. niger and N. hubbsi using the amh/Amh homolog of S. salar as an outgroup. CDS and proteins were predicted with the FGENESH+ suite (Solovyev et al., 2006), based on genomic sequences. These putative CDS and protein sequences were then aligned using MAFFT (version 7.450) (Katoh and Standley, 2013). Residue-wise confidence scores were computed with GUIDANCE 2 (Sela et al., 2015), and only well-aligned residues with confidence scores above 0.99 were retained. The resulting alignment file was used for model selection and tree inference with IQ-TREE (version 1.6.7) (Minh et al., 2020) with 1000 bootstraps and the 1000 SH-like approximate likelihood ratio test for robustness. To verify that the topology of the single amh homolog of D. pectoralis was not an artifact of long-branch attraction, we also constructed phylogenetic trees with only the first and second codons of the coding sequences (Lemey et al., 2009) as well as with complete protein sequences. For additional confirmation of the tree topology for the amh homologs, the same alignments were run in a Bayesian framework implemented in Phylobayes (version 4.1), (Lartillot et al., 2007; Lartillot and Philippe, 2006, Lartillot and Philippe, 2004) using the CAT-GTR model with default parameters, and two chains were run in parallel for approximately 1000 cycles until the average standard deviation of split frequencies remained ≤0.01.
Detection of signatures of selection on amha and amhby genes
Request a detailed protocolTo compare selection pressure between amha and amhby genes, dN/dS ratios were computed between each ortholog sequence and the amh sequence of U. pygmaea, which diverged from the other Esociformes prior to the amh duplication. Sequences from E. niger and N. hubbsi were excluded because their amhby orthologs were substantially shorter and could introduce bias in the analysis. The ratio of nonsynonymous to synonymous substitution (dN/dS = ) among amha and amhby sequences was estimated using PALM (version 4.7) (Yang, 2007) based on aligned full-length CDS, and the phylogenetic tree obtained with the CDS was used in the analysis. First, a relative rate test on amino acid substitution (Hughes, 1994) was performed between amha and amhby pairs with the amh of U. pygmaea used as an outgroup sequence. For each species with amh duplication, was calculated between the amha ortholog and amh of U. pygmaea and was compared to the calculated between the amhby ortholog and amh of U. pygmaea. A Wilcoxon test was used to compare the resulting values for the amh paralogs of these species. Then, several branch and site models (M0, M1a, M2a, M7, M8, free-ratio, and branch-site) implemented in the CODEML package were fitted to the data. The goodness of fit of these models was compared using the likelihood ratio test implemented in PALM.
Genome assembly
Request a detailed protocolRaw Illumina sequencing reads for D. pectoralis, E. masquinongy, E. niger, and N. hubbsi were assembled using the DISCOVAR de novo software with the following assembly parameters: MAX_MEM_GB = 256, MEMORY_CHECK = False, and NUM_THREADS = 16. Assembly metrics were calculated with the assemblathon_stats.pl script (Earl et al., 2011), and assembly completeness was assessed with BUSCO (version 3.0.2) (Simão et al., 2015) using the Actinopterygii gene set (Supplementary file 10).
To facilitate the comparison of sex loci in different populations of E. lucius, we improved the continuity of the previously published assembly of an E. lucius European male (GenBank assembly accession: GCA_007844535.1). This updated assembly was performed with Flye (version 2.6) (Kolmogorov et al., 2019) using standard parameters and genome-size set to 1.1 g to match theoretical expectations (Hardie and Hebert, 2004). Two rounds of polishing were performed with Racon (version 1.4.10) (Vaser et al., 2017) using default settings and the Nanopore reads aligned to the assembly with minimap2 (version 2.17) (Li, 2018) with the ‘map-ont’ preset. Then, three additional rounds of polishing were performed with Pilon (version 1.23) (Walker et al., 2014) using the ‘fix all’ setting and the Illumina reads aligned to the assembly with BWA mem (version 0.7.17) (Li and Durbin, 2009). Metrics for the resulting assembly were calculated with the assemblathon_stats.pl script (Earl et al., 2011). The completeness of the assembly was assessed with BUSCO (version 3.0.2) (Simão et al., 2015) using the Actinopterygii gene set (4584 genes) and the default gene model for Augustus. The resulting assembly was scaffolded with RaGOO (version 1.1) (Alonge et al., 2019) using the published female genome (GenBank accession: GCA_011004845.1) anchored to chromosomes as a reference.
Appendix 1
Data availability
All gene sequences, genomic, Pool-seq and RAD-Seq reads were deposited under the common project number PRJNA634624.
-
NCBI BioProjectID PRJNA634624. Sex determination in the Esociformes.
-
NCBI NucleotideID GCA_000721915.2. Esox lucius isolate CL-BC-CA-002, whole genome shotgun sequencing project.
References
-
Gapped BLAST and PSI-BLAST: a new generation of protein database search programsNucleic Acids Research 25:3389–3402.https://doi.org/10.1093/nar/25.17.3389
-
Environmental effects on fish sex determination and differentiationSexual Development 3:118–135.https://doi.org/10.1159/000223077
-
A duplicated, truncated amh gene is involved in male sex determination in an old world silversideG3: Genes, Genomes, Genetics 117:042697.https://doi.org/10.1534/g3.117.042697
-
Stacks: an analysis tool set for population genomicsMolecular Ecology 22:3124–3140.https://doi.org/10.1111/mec.12354
-
Processing of anti-mullerian hormone regulates receptor activation by a mechanism distinct from TGF-betaMolecular Endocrinology 24:2193–2206.https://doi.org/10.1210/me.2010-0273
-
Pfam: the protein families databaseNucleic Acids Research 42:D222–D230.https://doi.org/10.1093/nar/gkt1223
-
Novel sex chromosomes in 3 cichlid fishes from lake tanganyikaJournal of Heredity 109:489–500.https://doi.org/10.1093/jhered/esy003
-
ExPASy: the proteomics server for in-depth protein knowledge and analysisNucleic Acids Research 31:3784–3788.https://doi.org/10.1093/nar/gkg563
-
BookSex Determination and Sex Control in SalmonidaeIn: Fostier A, editors. Sex Control in Aquaculture. John Wiley & Sons. pp. 249–280.https://doi.org/10.1002/9781119127291.ch11
-
Genome-size evolution in fishesCanadian Journal of Fisheries and Aquatic Sciences 61:1636–1646.https://doi.org/10.1139/f04-106
-
Genotypic sex determination in teleosts: insights from the testis-determining amhy geneGeneral and Comparative Endocrinology 192:55–59.https://doi.org/10.1016/j.ygcen.2013.03.019
-
Predicting RAD-seq marker numbers across the eukaryotic tree of lifeGenome Biology and Evolution 7:3207–3225.https://doi.org/10.1093/gbe/evv210
-
The evolution of functionally novel proteins after gene duplicationProceedings. Biological Sciences 256:119–124.https://doi.org/10.1098/rspb.1994.0058
-
MAFFT multiple sequence alignment software version 7: improvements in performance and usabilityMolecular Biology and Evolution 30:772–780.https://doi.org/10.1093/molbev/mst010
-
Novel sex-determining genes in fish and sex chromosome evolution: novel Sex-Determining genes in fishDevelopmental Dynamics 242:339–353.https://doi.org/10.1002/dvdy.23927
-
Assembly of long, error-prone reads using repeat graphsNature Biotechnology 37:540–546.https://doi.org/10.1038/s41587-019-0072-8
-
TimeTree: a resource for timelines, timetrees, and divergence timesMolecular Biology and Evolution 34:1812–1819.https://doi.org/10.1093/molbev/msx116
-
A bayesian mixture model for across-site heterogeneities in the amino-acid replacement processMolecular Biology and Evolution 21:1095–1109.https://doi.org/10.1093/molbev/msh112
-
Computing bayes factors using thermodynamic integrationSystematic Biology 55:195–207.https://doi.org/10.1080/10635150500433722
-
The sequence alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
Minimap2: pairwise alignment for nucleotide sequencesBioinformatics 34:3094–3100.https://doi.org/10.1093/bioinformatics/bty191
-
Weird animal genomes and the evolution of vertebrate sex and sex chromosomesAnnual Review of Genetics 42:565–586.https://doi.org/10.1146/annurev.genet.42.110807.091714
-
IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic eraMolecular Biology and Evolution 37:1530–1534.https://doi.org/10.1093/molbev/msaa015
-
The dynamics and evolutionary potential of domain loss and emergenceMolecular Biology and Evolution 29:787–796.https://doi.org/10.1093/molbev/msr250
-
Sexual antagonism and the instability of environmental sex determinationNature Ecology & Evolution 2:343–351.https://doi.org/10.1038/s41559-017-0427-9
-
Turnover of sex chromosomes in Celebensis Group Medaka FishesG3: Genes, Genomes, Genetics 5:2685–2691.https://doi.org/10.1534/g3.115.021543
-
BookMechanisms of Environmental Sex Determination in Fish, Amphibians, and ReptilesIn: Navara K. J, editors. Choosing Sexes: Mechanisms and Adaptive Patterns of Sex Allocation in Vertebrates, Fascinating Life Sciences. Springer International Publishing. pp. 213–240.https://doi.org/10.1007/978-3-319-71271-0_10
-
BookPredicting Secretory Proteins with SignalP. InProtein Function PredictionHumana Press.
-
EST-based microsatellites for northern pike (Esox lucius) and cross-amplification across all Esox speciesConservation Genetics Resources 6:451–454.https://doi.org/10.1007/s12686-013-0123-2
-
BookEvolution of Sex Determining Genes in Fish Encyclopedia of ReproductionLaboratoire de Physiologie et Génomique des Poissons.
-
Sex determination: why so many ways of doing it?PLOS Biology 12:e1001899.https://doi.org/10.1371/journal.pbio.1001899
-
Genetic structure of muskellunge in the great lakes region and the effects of supplementation on genetic integrity of wild populationsJournal of Great Lakes Research 43:1141–1152.https://doi.org/10.1016/j.jglr.2017.09.005
-
Molecular and evolutionary dynamics of animal sex-chromosome turnoverNature Ecology & Evolution 3:1632–1641.https://doi.org/10.1038/s41559-019-1050-8
-
BookFreshwater Fishes of North America: Volume 2: Characidae to PoeciliidaeJohns Hopkins University Press.
-
PAML 4: phylogenetic analysis by maximum likelihoodMolecular Biology and Evolution 24:1586–1591.https://doi.org/10.1093/molbev/msm088
Article and author information
Author details
Funding
Agence Nationale de la Recherche (ANR-13-ISV7-0005)
- Yann Guiguen
Deutsche Forschungsgemeinschaft
- Manfred Schartl
Agence Nationale de la Recherche (ANR-10-INBS-09)
- Celine Roques
- Laurent Journot
Korea National Institute of Health (R01GM085318)
- John Postlethwait
ANR (SCHA 408/10–1)
- Yann Guiguen
- Manfred Schartl
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are grateful to the fish facility of INRAE LPGP for support in experimental installation and fish maintenance and to the GenoToul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo GenoToul) for providing help, computing, and storage resources. This work was supported by the Agence Nationale de la Recherche, the Deutsche Forschungsgemeinschaft (ANR/DFG, PhyloSex project, 2014–2016, SCHA 408/10–1, to YG and MS), and the U.S. National Institutes of Health (R01GM085318 to JHP). The MGX and Get-Plage core sequencing facilities were supported by France Genomique National infrastructure, funded as part of the ‘Investissement d’avenir’ program managed by the Agence Nationale pour la Recherche (contract ANR-10-INBS-09). We thank the MNHN curators for providing access to the collections of E. aquitanicus and E. cisalpinus and for the ONEMA. We thank GB Delmastro for help with specimen sampling. We also thank Mackenzie Garvey and Penny Swanson for assistance with dissection of N. hubbsi.
Copyright
© 2021, Pan et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,341
- views
-
- 306
- downloads
-
- 30
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Evolutionary Biology
- Genetics and Genomics
Evolutionary arms races can arise at the contact surfaces between host and viral proteins, producing dynamic spaces in which genetic variants are continually pursued. However, the sampling of genetic variation must be balanced with the need to maintain protein function. A striking case is given by protein kinase R (PKR), a member of the mammalian innate immune system. PKR detects viral replication within the host cell and halts protein synthesis to prevent viral replication by phosphorylating eIF2α, a component of the translation initiation machinery. PKR is targeted by many viral antagonists, including poxvirus pseudosubstrate antagonists that mimic the natural substrate, eIF2α, and inhibit PKR activity. Remarkably, PKR has several rapidly evolving residues at this interface, suggesting it is engaging in an evolutionary arms race, despite the surface’s critical role in phosphorylating eIF2α. To systematically explore the evolutionary opportunities available at this dynamic interface, we generated and characterized a library of 426 SNP-accessible nonsynonymous variants of human PKR for their ability to escape inhibition by the model pseudosubstrate inhibitor K3, encoded by the vaccinia virus gene K3L. We identified key sites in the PKR kinase domain that harbor K3-resistant variants, as well as critical sites where variation leads to loss of function. We find K3-resistant variants are readily available throughout the interface and are enriched at sites under positive selection. Moreover, variants beneficial against K3 were also beneficial against an enhanced variant of K3, indicating resilience to viral adaptation. Overall, we find that the eIF2α-binding surface of PKR is highly malleable, potentiating its evolutionary ability to combat viral inhibition.
-
- Ecology
- Evolutionary Biology
Seasonal polyphenism enables organisms to adapt to environmental challenges by increasing phenotypic diversity. Cacopsylla chinensis exhibits remarkable seasonal polyphenism, specifically in the form of summer-form and winter-form, which have distinct morphological phenotypes. Previous research has shown that low temperature and the temperature receptor CcTRPM regulate the transition from summer-form to winter-form in C. chinensis by impacting cuticle content and thickness. However, the underling neuroendocrine regulatory mechanism remains largely unknown. Bursicon, also known as the tanning hormone, is responsible for the hardening and darkening of the insect cuticle. In this study, we report for the first time on the novel function of Bursicon and its receptor in the transition from summer-form to winter-form in C. chinensis. Firstly, we identified CcBurs-α and CcBurs-β as two typical subunits of Bursicon in C. chinensis, which were regulated by low temperature (10 °C) and CcTRPM. Subsequently, CcBurs-α and CcBurs-β formed a heterodimer that mediated the transition from summer-form to winter-form by influencing the cuticle chitin contents and cuticle thickness. Furthermore, we demonstrated that CcBurs-R acts as the Bursicon receptor and plays a critical role in the up-stream signaling of the chitin biosynthesis pathway, regulating the transition from summer-form to winter-form. Finally, we discovered that miR-6012 directly targets CcBurs-R, contributing to the regulation of Bursicon signaling in the seasonal polyphenism of C. chinensis. In summary, these findings reveal the novel function of the neuroendocrine regulatory mechanism underlying seasonal polyphenism and provide critical insights into the insect Bursicon and its receptor.