Parallel Evolution of Sex-Linked Genes across XX/XY and ZZ/ZW Sex Chromosome Systems in the Frog Glandirana rugosa

Genetic sex-determination features male (XX/XY) or female heterogamety (ZZ/ZW). To identify similarities and differences in the molecular evolution of sex-linked genes between these systems, we directly compared the sex chromosome systems existing in the frog Glandirana rugosa. The heteromorphic X/Y and Z/W sex chromosomes were derived from chromosomes 7 (2n = 26). RNA-Seq, de novo assembly, and BLASTP analyses identified 766 sex-linked genes. These genes were classified into three different clusters (XW/YZ, XY/ZW, and XZ/YW) based on sequence identities between the chromosomes, probably reflecting each step of the sex chromosome evolutionary history. The nucleotide substitution per site was significantly higher in the Y- and Z-genes than in the X- and W- genes, indicating male-driven mutation. The ratio of nonsynonymous to synonymous nucleotide substitution rates was higher in the X- and W-genes than in the Y- and Z-genes, with a female bias. Allelic expression in gonad, brain, and muscle was significantly higher in the Y- and W-genes than in the X- and Z-genes, favoring heterogametic sex. The same set of sex-linked genes showed parallel evolution across the two distinct systems. In contrast, the unique genomic region of the sex chromosomes demonstrated a difference between the two systems, with even and extremely high expression ratios of W/Z and Y/X, respectively.


Introduction
In vertebrates, sex is determined genotypically or environmentally [1,2]. In the genotypic sex-determination system, a sex-determining gene on the sex chromosome triggers the primary formation of the testis or ovary. The heterogametic sex can be male (XY) or female (ZW). In homeothermic vertebrates such as mammals and birds, the heterogametic sex is male and female, respectively. Each system has been highly conserved within its own taxon for hundreds of millions of years [3,4]. Conversely, in poikilothermic vertebrates, the heterogametic sex is different between taxa, species, or even geographic populations within a species, and a transition between the two systems is possible [5][6][7][8]. In both systems, the Y or W chromosome dominates, or a single dose of the Z or X chromosome determines one sex [1]. Here, a basic question arises about the two systems. What are the differences in the evolutionary strategies of the sex chromosomes and sex-linked genes between the two distinct systems? In mammals and birds, the mutation rates of the Y and Z chromosomes are higher than those of the X and W chromosomes, which are male-driven mutations, and nonsynonymous to synonymous substitution rate ratios are higher in the genes of Y and W fication were performed using BLAST+ 2.2.31 (BLASTP, ≥90% identity, ≥100 aa, https: //blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download, accessed on 18 February 2021) [31]. Using the BLASTP analysis, duplicated gene sequences were removed and the longest orthologs per gene were extracted. XX, XY, ZZ, or ZW RNA-sequencing was mapped to the assembled XX or ZZ CDSs using BWA 0.7.17 (https: //bio-bwa.sourceforge.net/, accessed on 1 March 2021) [32]. The duplicate reads were removed using Picard 2.25.0 (http://broadinstitute.github.io/picard/, accessed on 1 March 2021). We used GATK4 v4.2.0.0 (https://gatk.broadinstitute.org/hc/en-us/articles/36 0036194592-Getting-started-with-GATK4, accessed on 1 March 2021) to perform variant calling [33]. The sequences of the X, Y, Z, and W chromosome-linked genes were constructed using the "consensus" command in bcftools [34]. BLASTP analysis of the gene sequences between G. rugosa and R. temporaria indicated that sex chromosome 7 of G. rugosa is orthologous to chromosome 9 (autosome) of R. temporaria ( Figure S1B). The analysis scheme is shown in Figure S1A.

Phylogenetic Analyses of Sex Chromosome-Linked Sox3 and tl Genes
We have previously performed phylogenetic analyses using sex chromosome-linked nucleotide sequences of sox3 and tl genes in the XY (Sekigahara and Hamamatsu cities), ZW (Kanazawa and Niigata cities), West, and East groups [25]. To confirm the quality of the sox3 and tl gene models in this study, we constructed phylogenetic trees based on the constructed gene sequences including the previous dataset [25] by using the neighbor joining (NJ) method ( Figure S4). The NJ trees of the two genes showed similar topologies comprising two major clades. One clade contained the genes from autosome 7 of the East-J group (Chiba city), X chromosomes of the XY group (Sekigahara, Hmamatsu, and Ichinomiya cities), and W chromosomes of the ZW group (Kanazawa, Niigata, and Suzu cities). The other clade contained the genes of autosome 7 of the West-J group (Hiroshima city), Y chromosomes of the XY group, and Z chromosomes of the ZW group. These results are consistent with the evolutionary scenario of the sex chromosomes: Y and Z chromosomes were derived from autosome 7 of the West-Japan group, whereas X and W chromosomes were derived from autosome 7 of the East-Japan group (G. reliquia) ( Figure 1A) [23,25]. It was then confirmed that the gene models constructed in this study have sufficient qualities for analysis on the molecular evolution of the X, Y, Z and W chromosome-linked genes in this species.

Expression Analyses
The analysis of the allele-specific expression between the X and Y chromosomelinked genes in the XY male or Z and W chromosome-linked genes in the ZW female was performed based on each sex chromosome-linked genes using ASE-TIGAR (http: //nagasakilab.csml.org/ase-tigar/, accessed on 17 June 2021) [43]. The statistical analysis was performed using Wilcoxon-Mann-Whitney tests.

Identification of X, Y, Z, and W Sex Chromosomes-Linked Genes
To identify the sex chromosome-linked coding sequence of G. rugosa, we performed RNA-sequencing, de novo assembly, and gene annotation (see "Section 2"). We then blasted the identified RNA sequences to the R. temporaria database and obtained 13,884 and 14,403 genes from the XX and ZZ RNA sequences, respectively. Of them, 12,572 genes were determined to be common orthologs between the XX and ZZ RNAs by BLASTP. Then, we mapped the XX and XY or ZZ and ZW RNA sequences to the 12,572 XX or ZZ genes, respectively. Finally, we constructed the X, Y, Z, and W sex chromosome-linked gene models, which comprised 766 sex chromosome-linked genes, 11,628 autosome-linked genes, and 178 unplaced genes ( Figure S1 and Table S1).
To confirm the sex linkage of the 766 sex chromosome-linked genes, concatenated phylogenetic trees of 11,628 autosomal chromosome-linked genes and the 766 genes were constructed by RaxML-NG. The trees displayed their evolutionary scenarios of autosomes and sex chromosomes, respectively ( Figures 1A and S2). The bootstrap values were 100 in the autosomal tree, but relatively lower (65%) in the sex chromosomal tree. The findings verified the Y-, X-, Z-, and W-genes as the sex chromosome-linked genes.

Three Clusters of Sex-Linked Genes
Based on the topology of the tree of each gene constructed by RaxML-NG, the 766 sexlinked genes were classified into three different clusters: XW/YZ (542 genes), XY/ZW (203 genes), and XZ/YW (21 genes) (Table S1). Then, the concatenated tree of each of the three clusters was constructed by RaxML-NG ( Figure 1B-D). The local bootstrap values in the three gene trees were 100, 100, and 80/93, respectively ( Figure 1B-D).
We then investigated the genomic distributions of the genes belonging to the three clusters on chromosome 9 of R. temporaria, which is orthologous to the sex chromosomes of G. rugosa. The genes of the three clusters were distributed evenly along the chromosomal axis, except for one region spanning from 160 to 184 M (out of 0-184 M), where 84.6% of the genes of the XY/ZW cluster were concentrated, while 10.8% and 4.6% of the genes of the XW/YZ and XZ/YW clusters, respectively, were concentrated ( Figure 2). The sequence identity between the Z-and W-genes of the XY/ZW cluster in the 160-184 M region was 99.87% and was significantly higher than 99.29% in the 0-160 M region (Tukey HSD, p < 0.001 in Figure 3C). On the other hand, the sequence identities between the X-and Y-genes were 99.22% and 99.47% in the 0-160 M and 160-184 M regions, respectively, which were not significantly different from each other (Turkey HSD, p = 0.1274, in Figure 3C). In contrast, the sequence identity between the X-and Y-genes of the XZ/YW cluster was 95.74% in the 160-184 M region and significantly lower than 99.32% in the 0-160 M region (Tukey HSD, p < 0.001 in Figure 3D).

Male-Biased Mutation
To elucidate the evolutionary rates of the X-, Y-, Z-, and W-genes, we performed a phylogenetic analysis of the 766 sex-linked genes of three clusters using the PAML free ratio model. We obtained the averaged values of nucleotide substitutions per codon (t) of the sex-linked genes in the branches from the nodes to tips of the X-, Y-, Z-, and W-genes.
The nucleotide substitution rates of the Y-and Z-genes were significantly higher than those of the X-and W-genes, respectively (Fisher exact test, p < 0.05, XY; p < 0.001, ZW) ( Figure 4).

Figure 2.
Composition of the sex-linked genes belonging to three different clusters along the chromosomal axis. The genes of the XW/YZ, XY/ZW, and XZ/YW clusters are indicated in blue, orange, and gray, respectively.

Male-Biased Mutation
To elucidate the evolutionary rates of the X-, Y-, Z-, and W-genes, we performed a phylogenetic analysis of the 766 sex-linked genes of three clusters using the PAML free ratio model. We obtained the averaged values of nucleotide substitutions per codon (t) of the sex-linked genes in the branches from the nodes to tips of the X-, Y-, Z-, and W-genes. The nucleotide substitution rates of the Y-and Z-genes were significantly higher than those of the X-and W-genes, respectively (Fisher exact test, p < 0.05, XY; p < 0.001, ZW) ( Figure 4).

Female-Biased dN/dS Ratio
To elucidate the strength of natural selection in the sex-linked genes, we calculated the number of synonymous and nonsynonymous nucleotide substitutions of the X-, Y-, Z-, and W-genes by PAML and then the number of nonsynonymous substitutions per nonsynonymous site (dN), number of synonymous substitutions per synonymous site (dS), and their ratios (dN/dS) of each gene (Figure 4) or concatenated sequences of each of the X-, Y-, Z-, and W-genes ( Table 1). The dN/dS X and dN/dS W ratios were higher than the dN/dS Y and dN/dS Z ratios, respectively, indicating female-biased dN/dS ratios (Table 1 and Figure 4). Subsequently, we performed Gene Ontology enrichment analysis in the four different categories of dN/dS ratios (dN/dS Y > dN/dS X , dN/dS Y < dN/dS X , dN/dS Z > dN/dS W , and dN/dS Z < dN/dS W ) in the sex chromosome linked genes. The categories did not differ significantly from each other. No selection pressure on gene function at the chromosome scale was observed during sex chromosome evolution in this species.
We also investigated the distribution of t, dN, dS, and dN/dS (X and Z values were subtracted from Y and W values, respectively) of the sex-linked genes along the chromosomal axis and compared the two systems. t, dS, and dN/dS were negatively correlated (R = −0.40, −0.51, and −0.50, respectively), whereas dN was slightly positively correlated with each other (R = 0.17) ( Figure 5).

Y-and W-Biased Expressions
To elucidate the allelic expression of the sex-linked genes in different organs of G. rugosa, we calculated the allele-specific expression of each of the three clusters between the X, Y, Z, and W chromosomes. Y-and W-biased expression was observed in the XW/YZ and XY/ZW clusters (p < 0.05 or p < 0.001), but not in the XZ/YW cluster, except in the ZW brain (p < 0.01) ( Figure 6 and Table S2).
Next, we investigated the distribution of the Y/X and W/Z expression ratios along the chromosomal axis. The two expression ratios were positively correlated with each other at the whole region from 0-184 M (R = 0.22). Particularly unique was that at the position 160-184 M, the Y/X and W/Z ratios were negatively correlated with each other (R = −0.42) in contrast to that in the 0-160 M region (R = 0.36), where the Y/X ratios were very high, while the W/Z ratios were almost even (Figures 7 and S3).

Y-and W-Biased Expressions
To elucidate the allelic expression of the sex-linked genes in different orga rugosa, we calculated the allele-specific expression of each of the three clusters the X, Y, Z, and W chromosomes. Y-and W-biased expression was observed in the and XY/ZW clusters (p < 0.05 or p < 0.001), but not in the XZ/YW cluster, except in brain (p < 0.01) ( Figure 6 and Table S2).
Next, we investigated the distribution of the Y/X and W/Z expression ratios a chromosomal axis. The two expression ratios were positively correlated with ea at the whole region from 0-184 M (R = 0.22). Particularly unique was that at the

Three Evolutionary Strata of the Sex Chromosomes in G. rugosa
Based on the topologies of the phylogenetic trees of the sex-linked genes, we fied three clusters: XW/YZ (542 genes), XY/ZW (203 genes), and XZ/YW (21 gene ures 1 and 2). The first cluster, XW/YZ, may have been derived from the autosoma mosome 7 of the ancestral populations before their hybridization. Thus, this stratu be the oldest. The XY and ZW sex chromosome systems in G. rugosa share a phylo origin of past hybridization between the populations of West-Japan and East-Jap reliquia) [5,23,24]. Evidence can be seen from the phylogenetic origins of t

Three Evolutionary Strata of the Sex Chromosomes in G. rugosa
Based on the topologies of the phylogenetic trees of the sex-linked genes, we identified three clusters: XW/YZ (542 genes), XY/ZW (203 genes), and XZ/YW (21 genes) (Figures 1 and 2). The first cluster, XW/YZ, may have been derived from the autosomal chromosome 7 of the ancestral populations before their hybridization. Thus, this stratum may be the oldest. The XY and ZW sex chromosome systems in G. rugosa share a phylogenetic origin of past hybridization between the populations of West-Japan and East-Japan (G. reliquia) [5,23,24]. Evidence can be seen from the phylogenetic origins of the sex

Three Evolutionary Strata of the Sex Chromosomes in G. rugosa
Based on the topologies of the phylogenetic trees of the sex-linked genes, we identified three clusters: XW/YZ (542 genes), XY/ZW (203 genes), and XZ/YW (21 genes) (Figures 1 and 2). The first cluster, XW/YZ, may have been derived from the autosomal chromosome 7 of the ancestral populations before their hybridization. Thus, this stratum may be the oldest. The XY and ZW sex chromosome systems in G. rugosa share a phylogenetic origin of past hybridization between the populations of West-Japan and East-Japan (G. reliquia) [5,23,24]. Evidence can be seen from the phylogenetic origins of the sex chromosomes based on the sex-linked gene sequences [23,25]. In this study, we confirmed the dual origin of sex chromosomes ( Figure 1B). Based on the gene trees that included orthologous genes of two ancestral populations, it is evident that the Z/Y chromosomes derived from the orthologous autosome 7 of the West-Japan group, while the W/X chromosomes derived from autosome 7 of the East -Japan group (G. reliquia) ( Figure S4).
The second cluster of XY/ZW may have had two different origins. One is a protosex chromosomal origin; the cluster was generated after hybridization between the two ancestral populations and just before or at the sex chromosomal establishment ( Figure S5). The two different autosomes 7 might have recombined with each other in the male and female meioses of the hybrids. After the two sex chromosome systems were established and separated into different populations, they started to accumulate independent nucleotide substitutions in each of the two systems. The other is a pseudoautosomal region (PAR) origin generated in the ZW system only ( Figure S5). The sequence identities between the Zand W-genes in the 160-184 M region were significantly higher than those from the other 0-160 M region. Allelic expression of the genes from the 160-184 M region also did not show any W-bias. Thus, the 160-184 M region may be a PAR of the ZW sex chromosomes. Morphologically, the terminal regions of the Z chromosome short arm and W chromosome long arm are estimated to be PARs based on the distribution of chiasmata during female meiosis and the shared signals in the male and female karyotypes by comparative genomic hybridization analysis [44,45]. In contrast, PAR was not identified in the XY system; under a microscope, the X and Y chromosomes were observed to be paired by end-to-end formation with no chiasmata during male meiosis [25]. In addition, in highly evolved frogs such as true frogs including G. rugosa, the bivalent chromosomes in male meiosis form a ringshaped bivalent with no internal chiasma. On the other hand, those in female meiosis form chiasmata along the chromosome axis, except around the centromeric region [46][47][48]. Thus, the PAR of the X and Y chromosomes, if any, may be an extremely tiny terminal tip region and may differ markedly in size from the PAR of the ZW system.
The third cluster of XZ/YW may have been generated by hybridization between the two systems ( Figure S6); this stratum may be the youngest. Our previous study showed that XY populations emigrated into the ZW populations after the two systems were established [49], and that the X and Z chromosomes or Y and W chromosomes were recombined with each other because the WY and XZ frogs can be generated by crossing between the frogs from the two systems [50].
The three strata in G. rugosa were structurally different from those observed in the sex chromosomes of other vertebrates, in which one stratum is built by chromosomal rearrangement such as inversion, and the gene members of each cluster are arranged in a sequence [51][52][53]. In contrast, in G. rugosa, the strata were dispersed entirely on the chromosomes, except for the PAR of the ZW chromosomes, and may have been built on the process of sex chromosome evolution from ancestral autosomes through proto-sex chromosomes to the third phase of inter-population hybridization after sex chromosome establishment [49]. The second and third strata may have formed through recombination between homologous chromosomes.

Male-Driven Mutation
Male-biased mutations in the Y-and Z-borne genes have been observed in many vertebrates and invertebrates [12,[54][55][56]. DNA replication errors during germline cell division are a major source of mutations that are transmitted to the next generation. Because male germ cells proliferate more frequently than female germ cells, nucleotide substitution rates tend to be higher in males than in females [57]. Y chromosomes are always present in males, and Z chromosomes are carried during two-thirds of the lives of generations by males [54]. Therefore, Y and Z chromosome-linked genes evolve faster than X and W chromosome-linked genes. In our previous study, a male-driven mutation was shown to work in the XY and ZW systems of G. rugosa based on several sex-linked genes [25]. The comprehensive analysis in this study strongly supports the previous results by showing male-biased substitution rates in the XY and ZW systems.

Constraints by Negative Selection on Sex-Linked Genes
The dN/dS ratio is a measure of the strength of natural selection acting on proteincoding genes [41]. The Y and W chromosomal genes show higher dN/dS ratios than the X and Z chromosomal genes in mammals [9,10], lizards [58], Drosophila [11,59], Silene [13,14], and birds [15]. Many of these are considered to be due to deleterious mutations occurring in the non-recombining regions of the Y and W chromosomes, leading to degeneration, or inversely, in some specialized genes under positive selection [15,60]. In G. rugosa, XY and ZW sex chromosomes are heteromorphic in males and females, respectively [18]. Two inversions, one in Z/Y chromosomes and the other in W/X chromosomes, each from primordial autosome 7, created non-recombining regions on the XY and ZW sex chromosomes [5] and the pseudoautosomal regions of the ZW chromosomes were restricted to their terminals [25,45,61]. Evidently, degeneration of the Y and W chromosomes is progressing because artificially constructed YY and WW embryos die of edemata at early developmental stages due to lethal, degenerated genes responsible for development on the Y and W chromosomes [49,62]. Therefore, we expected that the dN/dS ratios would be biased toward the Y-and W-borne genes in this species, as in mammals and birds. However, the dN/dS in all four sex chromosome-linked genes was small and under 1, indicating that purifying selection acts to conserve the functions of the sex-linked genes. In addition, dN/dS was higher in the W-genes (0.4469) than in the Z-genes (0.3473) in the ZW system, and was higher in the X-genes (0.3246) than in the Y-genes (0.2633) in the XY system. Provean scores, which can estimate the type of selection acting on a gene and are lower in deleterious mutations [42], were totally higher in the Z-and X-genes than in the W-and Y-genes ( Figure S7). These scores predict less deleterious substitutions or positive selection in Z-and X-genes. A similar case was observed in the stickleback fish Pungitius pungitius, in which the sex chromosomes were evolutionarily young and the dN/dS and Provean scores were larger in the X chromosome than in the Y chromosome [63].

Y-and W-Biased Expression
The allele expression was evidently higher in the Y-and W-genes than in the X-and Zborne genes, respectively, in the two clusters of XW/YZ and XY/ZW of G. rugosa. This result is in sharp contrast to other animals where Y-or W-borne genes showed lower expression than their homologues due to the accumulation of deleterious mutations [64][65][66]. Negative correlations between dN/dS ratios and expression levels have been observed [67,68]. We also found that the dN/dS ratios of the W-and Y-genes in G. rugosa were negatively correlated between the two systems along the chromosomal axis, while the Y/X and W/Z expression ratios were positively correlated (Figures 5 and 7). These correlation patterns suggest that dN/dS is not related to the allelic expression of sex-linked genes in either system (R = 0.16 in the XY system and R = 0.05 in the ZW system). We are likely to observe sex-linked gene expression patterns favoring heterogametic sex at the primary stage of young sex chromosome evolution. Likewise, the young sex chromosomes of stickleback fish (Pungitius pubgitius) showed higher Y/X expression ratios with lower rates of dN/dS in Y-borne genes [63].
Why is the expression higher in Y-and W-borne genes favoring heterogametic sex? As the data shown above are average values, we investigated the relationships between the Provean scores and allelic expression ratios for every gene from the three clusters. Y-genes with lower Provean scores than X-genes were much higher in expression than the X-genes, while the W-genes with higher Provean scores than Z-genes were much higher in expression than the Z-genes, except in the gonads (Figures S8-S10). These results suggest that less deleterious mutations or positive selection to W-genes are involved in the upregulation of the W-genes, but reason for the Y-gene upregulation may differ and is still unclear.

Unique Terminal Region of the Sex Chromosomes
The Y/X and W/Z expression ratios were positively correlated with each other along the chromosomal axis in the 0-160 M region. In contrast, at the 160-184 M terminal region, the two expression ratios were rather negatively correlated (Figure 7). This may be the PAR of the ZW sex chromosomes, as described above, and a tiny terminal PAR in the XY chromosomes. Although the W/Z expression ratios were nearly even, the Y/X expression ratios were remarkably high in this region. The 24 M terminal region may be a critically different region between the XY and ZW sex chromosomes in this species as well as a good genomic region for evolving a male-determining gene with a much higher Y/X expression rate in the XY system. Next, we searched for candidate genes to determine males with higher Y/X ratios among the regions (Table 2), where one candidate was identified, a steroidogenic enzyme, 17β-hydroxysteroid dehydrogenase 1 (Hsd17B8), which catalyzes the production of estradiol for ovarian differentiation. The W/Z expression ratios of the gene in the three tissues were approximately 1, while the Y/X ratios were 2.58-5.16 (Table 2). HSD18B1, another gene of this family, has been identified as a female-determining gene in Seriola fishes [69], and HSD18B8 itself has been isolated as a differentially methylated gene expressed in the gonads of turtles at the temperature sex-determining stages [70]. A truncated form of Hsd17B8 on the Y chromosome might be involved in the testis determination in G. rugosa. Unexpectedly, at almost the terminal tip of the region, we identified another gene, dachshund family transcription factor 2 (Dach2), of which the W/Z expression ratios were extremely high, 50.4-707, while the Y/X ratios were 0.93-2.83. The orthologous gene in humans is reported to be a candidate for premature ovarian failure [71], and thus could be a candidate for determining females in the ZW system of G. rugosa. On the other hand, from the other genomic regions showing higher expression ratios of both Y/X and W/Z along the chromosomal axis, two genes, negative elongation factor complex member E (NELFE) located at the 31 M region and E3 SUMO-protein ligase ZBED1-like located in the 109 M region, showed very high ratios of both Y/X and W/Z. NELFE is expressed in the gonads of X. tropicalis [72], and thus may have the potential to regulate both male and female determination in the two systems. Consequently, to examine the sex-determining functions of the genes listed above, it will be necessary to examine their expression in undifferentiated gonads in genetic male and female tadpoles and then perform a functional analysis.

Conclusions
The frog G. rugosa distributed on the Japanese islands is suitable for use in investigating the similarities and differences in the molecular evolution of X, Y, Z, and W chromosome-linked genes derived from the same homologous chromosomes at the early stage of sex chromosome differentiation. We identified three clusters of sex-linked genes that illustrate the evolutionary strata in the history of sex chromosomes. Molecularly, we confirmed male-biased mutation, female-biased dN/dS ratio (all < 1), and Y-and W-gene-biased expression in the sex chromosome-linked genes, showing parallel evolution across the two distinct systems. Importantly, we identified a unique genomic region at the terminal part of the chromosomes, which may be the PAR of the ZW sex chromosomes. If present, the PAR in the XY sex chromosomes was very small and the expression of the Y-genes was much higher than that of the X-genes. This region may represent the genomic and expression differences between the two systems. In future studies, the analyses of X, Y, Z, and W chromosome-linked genes will be extended to orthologous genes that are still autosomal in the two ancestral populations [73]. The results could reveal the dynamics of molecular evolution from autosomes to sex chromosomes and vice versa, associated with sex chromosome turnover.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes14020257/s1. Figure S1: RNA-sequencing, de novo assembly, and gene annotation in Glandirana rugosa. (A) Schematic diagram of the analysis scheme in this study. XY group (XX female and XY male) and ZW group (ZZ male and ZW female) of G. rugosa collected from Ichinomiya city, Aichi Prefecture (XY group) and Suzu city, Ishikawa Prefecture (ZW group). The RNAs were used for sequencing, de novo assembly, and gene annotation. (B) Chromosome comparison between G. rugosa and Rana temporaria. The karyotypes of G. rugosa and R. temporaria are composed of common 13 chromosome pairs (2n = 26) [74]. The subchromosomal locations of 12 genes in G. rugosa were taken from previous studies [24,[75][76][77]. The gene information in R. temporaria was obtained from the WGS Project (GCF_905171775.1, aRanTem1). Sex chr. denotes sex chromosome. Figure S2: Phylogenetic analyses of sex-linked and autosomal genes. (A,B) Concatenated phylogenetic trees by rAxML-NG using 11,628 autosomal chromosome-linked genes (A) and 766 sex chromosomelinked genes (B). The best-fit models of nucleotide substitution were selected by Modeltest-NG. Numbers at each node denote the bootstrap percentage values based on 1000 replicates. Rana temporaria was used as an outgroup. Figure S3: Comparison of the sex-linked gene expression between the sex chromosomes located at the 160-184 M regions. Statistical analysis of the expression of the sex-linked genes between the sex chromosomes were performed using Wilcoxon-Mann-Whitney tests (paired). Figure S4: Phylogenetic trees of the sex chromosome-linked sox3 (A) and tl (B) genes in the XY and ZW groups. The corresponding sequences of the G. rugosa West-and East-groups, and of Pelophylax nigromaculatus as an outgroup, were also used. The neighbor joining phylogenetic trees were constructed by MEGA11 using the maximum composite likelihood method. Numbers at each node denote the bootstrap percentage values based on 1000 replicates. The X-and W-linked genes are boxed in red, the Y-and Z-linked genes and genes on homologous autosome 7 of West-Japan are in blue, and those on homologous autosome 7 of East-Japan (G. reliquia) are in orange. The X/W, Y/Z chromosomes and autosome 7 of West-Japan, and autosome 7 of East-Japan (G. reliquia) are shown in red, blue, and orange, respectively, and are put on the side of the gene tree. Chiba East 7, East-Japan group (G. reliquia) (Chiba city, chromosome 7); Hamamatsu X, XY group (Hamamatsu city, X chromosome); Hamamatsu Y, XY group (Hamamatsu city, Y chromosome); Hiroshima west 7, West group (Hiroshima city, chromosome 7); Ichinomiya X, XY group (Ichinomiya city, X chromosome); Ichinomiya Y, XY group (Ichinomiya city, Y chromosome); Kanazawa Z, ZW group (Kanazawa city, Z chromosome); Kanazawa W, ZW group (Kanazawa city, W chromosome); Niigata Z, ZW group (Niigata city, Z chromosome); Niigata W, ZW group (Niigata city, W chromosome); Sekigahara X, XY group (Sekigahara city, X chromosome); Sekigahara Y, XY group (Sekigahara city, Y chromosome); Suzu Z, ZW group (Suzu city, Z chromosome); Suzu W, ZW group (Suzu city, W chromosome). Figure S5: Schematic representation showing the birth of sex-linked genes belonging to the XY/ZW cluster. Autosomes 7 W and 7 E may have been recombined, as indicated in grey, after hybridization between the ancestral type-populations of West-Japan and East-Japan (G. reliquia). The genes located on the recombined regions accumulated independently of nucleotide substitutions after the establishment of the sex chromosomes. The PAR of the terminal regions of the Z-and W-chromosomes are shown in pale purple. Figure S6: Schematic representation showing the birth of sex-linked genes belonging to the XZ/YW cluster. XY populations have emigrated into the ZW populations in the past [49]. It is hypothesized that recombination between the X and Z in the XZ hybrids or Y-and W-chromosomes in WY hybrids occurred to produce the sex-linked genes belonging to the XZ/YW cluster. Figure S7: Provean scores in the sex-linked genes. Average value of Provean score was higher in the X-and Z-linked genes than in the Y-and W-linked genes. Statistical analysis was performed using Wilcoxon-Mann-Whitney tests (paired). Figure S8: Expression of sex-linked genes from the XW/YZ cluster and Provean scores. Expression of the Y-genes and W-genes with lower and higher Provean scores, respectively, was higher than those of the X-genes and Z-genes except in the XY gonad. Statistical analysis on the expression of the sex-linked genes between sex chromosomes was performed using Wilcoxon-Mann-Whitney tests (paired). Figure S9: Expression of sex-linked genes from the XY/ZW cluster and Provean scores. Expression of the Y-genes and W-genes with lower and higher Provean scores, respectively, was higher than those of the X-genes and Z-genes except in the ZW gonad. Statistical analysis on the expression of the sex-linked genes between sex chromosomes was performed using Wilcoxon-Mann-Whitney tests (paired). Figure S10: Expression of sex-linked genes from the XZ/YW cluster and Provean scores. Expression of the sex-linked genes was not related to the Provean scores. Statistical analysis on the expression of sex-linked genes between the sex chromosomes was performed using the Wilcoxon-Mann-Whitney tests (paired). Table S1: List of sex-linked genes from three different clusters. Table S2: Allelic expression of the sex-linked genes based on the data of RNA-Seq analysis in three tissues from the XY male and XX female in the XY system and ZZ male and ZW females in the ZW system. Informed Consent Statement: Not applicable.

Data Availability Statement:
The assembly of the Rana temporaria genomic data is available at https://github.com/DanJeffries/Rana-temporaria-genome-assembly/wiki, accessed on 24 May 2021 [6]. The gene data of Glandirana rugosa other than the Supplementary Materials presented in this study are available upon request from the corresponding authors.