Gene traffic mediated by transposable elements shaped the dynamic evolution of ancient sex chromosomes of varanid lizard

19 Lizards exhibit rapid turnovers and a much greater diversity of sex determination mechanisms 20 compared to birds and eutherians. This makes the conserved ZW sex chromosomes of 21 anguimorph lizards originated over 115 million years ago a seeming exception. We recently 22 discovered in an anguimorph lizard Varanus acanthurus (Vac) that its entire W chromosome 23 (chrW), but not chrZ is homologous to part of the chr2 by cytogenetic mapping, suggesting a 24 complex history of its sex chromosome evolution yet to be elucidated. To address this, we 25 assembled a chromosome-level genome, and provided evidence that the Vac sex chromosome 26 pair had undergone at least three times of recombination suppression, forming a similar pattern 27 of ‘evolutionary strata’ to that of birds or mammals. We identified the putative sex-determining 28 genes in the oldest evolutionary stratum that had first lost recombination. Comparison to other 29 lizard genomes dated the stepwise propagation of specific retrotransposon subfamilies shared 30 by chrW and chr2 to the varanid ancestor. These retrotransposons are also enriched near the 31 duplicated genes shared by the two chromosomes and probably mediated the recruitment of 32 many autosomal genes that rejuvenated the degenerating chrW, including members of a large 33 vomeronasal chemosensory receptor gene family V2R. Our results challenge the canonical 34 model of sex chromosome evolution, and suggest that the W or Y chromosome as a refugium of 35 repetitive elements, may recurrently recruit short-lived functional genes responsible for sexual 36 dimorphisms during its long-term course of degeneration. 37 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 2

38 Lizards harbor a rich diversity of sex determination (SD) mechanisms that is unparalleled by 39 other land vertebrates (Ezaz et al. 2009; Mezzasalma et al. 2021), providing great opportunities 40 for answering fundamental questions like why and how these different mechanisms emerge and 41 transit to each other, and how do sex chromosomes subsequently evolve. Such a diversity has 42 long been recognized by early systematic cytogenetic reports of male heterogametic (male XY, 43 female XX, like mammals), female heterogametic (male ZZ, female ZW, like birds), multiple 44 (e.g., X1X2Y in Calyptommatus lizard due to autosome fusion with the Y chromosome 45 (Yonenaga-Yassuda et al. 2005)) sex chromosomes in different species (Bull 1980), and 46 temperature sex determination mechanism (TSD) that co-exists with one of the genetic sex 47 determination (GSD) mechanisms in some species (e.g., spotted snow skink (Hill et al. 2018)). 48 Another aspect of the diversity is the diversified rate of SD turnovers between different lizard 49 lineages. A classic RAD-seq study identified between 17 to 25 SD transitions from 12 50 representative gecko species (Gamble et al. 2015). By contrast, a broad qPCR study 51 characterized an extraordinarily conserved ZW sex chromosome pair shared by almost all 52 studied Anguimorpha (Gila monster, beaded lizards, alligator lizards and varanids) species 53 except for slowworm (Rovatsos et al. 2019), and the very recently reported crocodile lizards 54 (Pinto et al. 2024). The latter strong evolutionary stasis seems to support the ‘evolutionary trap’ 55 hypothesis (Bull and Charnov 1985; Pokorná and Kratochvíl 2009) that heteromorphic sex 56 chromosomes resulting from recombination suppression precludes transitions into other GSD 57 systems, as suggested by sex chromosomes of birds and mammals. 58 However, the fine and complex course of lizard sex chromosome evolution can be 59 concealed due to the limited coverage and resolution of cytogenetic/qPCR/Illumina based 60 methods, particularly the absence of abundant sex-specific chrY or chrW sequences. Although 61 a recent burst of chromosome level lizard genomes produced by long-reads has characterized 62 the euchromatic chrX or chrZ (Westfall et al. 2021; Geneva et al. 2022; Koochekian et al. 2022; 63 Davalos-Dehullu et al. 2023; Leitão et al. 2023; Webster et al. 2024), little is known about their 64 often-heterochromatic homologs chrY or chrW that likely encompass the upstream SD genes. 65 We previously identified the candidate SD genes of an Anguimorpha species Varanus 66 acanthurus (Vac) using an Illumina-based scaffold-level genome (Zhu et al. 2022). Based on 67 this draft genome, we discovered a curious homology using fluorescence in situ hybridization 68 between the entire chrW, but not chrZ, and the tip of autosome chr2 (Figure 1a, termed chr1 in 69 the previous studies (Kinga and King 1975; Johnson et al. 2016; Iannucci et al. 2019), 70 Supplementary Note) (Dobry et al. 2025), using random probes specifically targeting the chrW. 71 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 3 This new finding questioned the reported simple history of deep evolutionarily static varanid sex 72 chromosomes that originated over 115 million years ago (Rovatsos et al. 2019). 73 Three hypothetical scenarios that are not mutually exclusive may explain the cytogenetic 74 homology: first, the chrW has acquired duplicated sequence segments from chr2 after it stopped 75 recombination with the chrZ. Second, the chr2 may have acquired duplicated sequences 76 specifically from the chrW. And the last, certain repetitive sequences have accumulated at both 77 the chr2 tip and the chrW but nowhere else in the genome, possibly because both regions have 78 low or no recombination. Any of these scenarios may result in gene movements between sex 79 chromosomes and autosomes that undergo very different evolutionary and molecular processes. 80 Previous studies reported gene movements between the chrX and autosomes (termed ‘gene 81 traffic’) in mammals (Emerson et al. 2004) and Drosophila (Betrán et al. 2002; Vibranovski et al. 82 2009). An excess of genes that move out of the chrX, producing a paucity of male-biased X-83 linked genes relative to autosomes can be the result of sexual antagonistic selection (Wu and 84 Xu 2003) that predicts the chrX will undergo demasculinization of gene expression (Sturgill et al. 85 2007) because it is more often inherited in the female than in male. It can be also explained by 86 other molecular processes including the meiotic X-inactivation (Turner 2007) , and dosage 87 compensation (Bachtrog et al. 2010) that specifically affect the sex chromosomes (Vicoso and 88 Charlesworth 2006). Individual cases of either direction of gene movement, but not a scale as 89 large as what we observed in Vac (Figure 1a), have also been reported on the human chrY 90 (Hughes et al. 2015; Xu et al. 2020) , the avian chrW (Xu et al. 2020), and the UV sex 91 determination system of kelp (Lipinska et al. 2017). Therefore, the Vac sex chromosomes 92 provide a model for elucidating the molecular and evolutionary mechanisms of large sequence 93 exchanges between sex chromosomes and autosomes. To test the three mentioned 94 hypotheses, in this work we produced a high-quality chromosome-level genome of a female Vac 95 individual, and inferred the evolutionary history of its sex chromosomes, particularly the origin of 96 the large sequence homology between the chrW and chr2. 97 98

99 A high-quality genome of Varanus acanthurus 100 About 100-fold PacBio long read sequences were produced in order to acquire nearly complete 101 genome and sex chromosome sequences. We anchored 97.3% of the genome (with an 102 assembled size of 1.55 Gb vs. the estimated size of 1.5 Gb (Zhu et al. 2022)) into 19 autosomes 103 and a chrZ by Hi-C chromatin contact data (Figure 1a, Supplementary Table S1). The high 104 quality of the genome can be evidenced by the large contig N50 length (114 Mb), high BUSCO 105 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 4 value (Figure 1b, 97.8% complete), and the annotated centromeric regions that are specifically 106 enriched for L1 elements relative to the rest of the chromosome region (Figure 1c, 107 Supplementary Fig. 1). Similar to other reptiles including birds, the 19 autosomes include 8 108 macrochromosomes, and 11 microchromosomes of smaller sizes (< 30Mb) but much higher 109 gene density (42 genes per Mb vs. 21 of macrochromosomes). 110 We identified the chrZ of 11.8 Mb long, with an expected 2-fold higher genomic read 111 depth in males than in females across most of the region (the sexually differentiated region, 112 SDR) lacking homologous recombination in females. The chrZ is homologous to the 113 microchromosome chr28 of chicken (Figure 1d), consistent with previous reports (Rovatsos et 114 al. 2019; Zhu et al. 2022). One chromosome end of pseudoautosomal region (PAR) of only 0.57 115 Mb long shows an equal read depth between sexes and likely retains the recombination. We 116 also identified a total of 418 scaffolds totaling 15 Mb, which can be mapped by female but not 117 male reads, as expected for chrW sequences (See Methods). Besides the annotated 118 centromeric regions, one chromosome end of chr2q and chr4, and the W-linked sequences 119 have a significantly (P < 0.05, Wilcoxon test) higher overall repeat content, particularly in certain 120 subfamilies of transposable elements (TEs, e.g., L1, Gypsy, Figure 1e, Supplementary Fig. 2), 121 than the rest genome (Figure 1c). Although the chrW is depleted (P < 0.05, Wilcoxon test) for 122 some families of DNA transposons and short interspersed nuclear elements (SINEs). This 123 supports the idea that sex-specific chromosomes act as a refugium for recently active TEs 124 (Peona et al. 2021) because they cannot be effectively purged by natural selection under a non-125 recombining environment (Charlesworth and Charlesworth 2000). 126 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 5 127 Figure 1 A new chromosome level genome and sex chromosomes of Varanus acanthurus 128 a) Hi-C contact map of the Varanus acanthurus (Vac) genome. The blue square indicates an 129 intact chromosome with preferential interactions (red dots) within rather than between each 130 chromosome. We show the karyotype of Vac below, and the DNA-FISH results using the probes 131 designed against the sequences of chrW (pink) and chrZ (green) (Dobry et al. 2025). b) The 132 snail plot shows the major statistics of genome assembly including the number and length of the 133 scaffolds, the length of the longest scaffold, the scaffold N50 and N90 values and the BUSCO 134 score based on the vertebrate protein (vertebrate_odb10) dataset. c) The bubble plot shows the 135 identification of sex chromosomes, with the bubble size scaled to the scaffold size. The chrW-136 linked scaffolds show a male-to-female ratio of mappable sites close to zero, whereas this ratio 137 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 6 is around 1 for chrZ and autosomes. In addition, chrZ shows nearly a two-fold male-to-female 138 read coverage ratio compared to the autosomes, while the ratio for the chrW remains close to 139 zero. The circos plot shows from the outer track to inner ones: the density of overall annotated 140 repetitive elements; DNA read coverage of male (blue) and female (red) individuals; density of 141 four classes of repetitive elements; and the GC content. All values are calculated in 100 kb 142 windows. d) Homology between the sex chromosomes of Vac and autosomes of chicken, with 143 each line representing one orthologous gene pair between the two species. e) The heatmap 144 shows the row-scaled TE content divided by class and family across chromosomes 145 (Macro/micro refers to mean values calculated from macro/micro chromosomes). 146 147 Varanid sex chromosomes have undergone three recombination suppression events 148 A shared feature of amniote sex chromosomes, regardless XY or ZW systems, is that they 149 usually underwent stepwise complete suppression of recombination, producing a punctuated 150 pattern of pairwise sequence divergence levels between adjacent sex-linked regions termed 151 ‘evolutionary strata’ (Vicoso et al. 2013; Cortez et al. 2014; Zhou et al. 2014). By change-point 152 analyses (see Methods), we demarcated the Z-linked SDR into three evolutionary strata (from 153 the old to young stratum termed as S0,1,2, Figure 2a), by their significantly (P < 0.05, Wilcoxon 154 test) and consistently sharp differences of ZW pairwise sequence identity (% of identical bases 155 per 100kb Z-linked sequences) and Z-linked male heterozygosity (SNP density calculated by 156 male reads aligning to the chrZ) levels compared to the adjacent strata. A significantly different 157 level of Z-linked heterozygosity between strata can be explained by the different time span that 158 each stratum region has experienced with reduction of effective population size due to that of 159 recombination (Charlesworth et al. 1987). This, to our knowledge, provides the first 160 demonstrated case of evolutionary strata in lizards. 161 Due to complete lack of recombination, only 10% and 12% of the S0 and S1 W-linked 162 genes (termed ‘gametologs’) are inferred to retain functions under our strict criterion (Materials 163 and Methods). The rest W-linked gametologs have either become completely deleted 164 compared to their Z-linked counterparts, or disrupted in the open reading frames (ORFs), or 165 silenced in expression in all examined female tissues (Figure 2b). This is a conserved estimate 166 of putatively functional genes on the chrW, as some ‘silenced’ genes are likely to be expressed 167 in other female tissues that are not included in this study. If we count those ‘silenced’ intact 168 ORFs, the putative functional W gametologs of Vac range between 21 to 30 genes, a 169 comparable number to that of human chrY. This translates to an average gene loss rate of 170 between 2.9 to 4.54 (0.53% to 0.87% of the) genes per million years on the chrW, a comparable 171 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 7 rate to that of human chrY (0.5% per million years (Krasovec et al. 2018)). Interestingly, the 172 resultant hemizygous Z-linked gametologs within the older S0, but not S1 and S2 show a 173 significantly faster evolution rate (Wilcoxon test, P < 0.05), measured by pairwise 174 nonsynonymous vs. synonymous substitution rate (dN/dS ratio) to the autosomal orthologs in 175 the desert horned lizard (Phrynosoma platyrhinos, Ppl) (Koochekian et al. 2022), than 176 autosomal genes, indicating a similar faster-Z effect (Charlesworth et al. 1987; Mank et al. 2010177 reported in birds (Wang et al. 2014) and some other female heterogametic sex systems (Vicoso 178 et al. 2013; Sackton et al. 2014) (Figure 2c). 179 We previously inferred the candidate sex determining gene of Vac to be the Z-linked 180 Amh, which has a high expression in the testis and male brain, but almost no expression in the 181 ovary and female brain (Zhu et al. 2022). Here we confirmed that Amh is located in the S0 182 region of chrZ, and we have not found its homolog on the chrW (Figure 2d). This is consistent 183 with the evolutionary hypothesis that recombination was first lost at the region encompassing 184 the sex determining gene, as demonstrated in birds and mammals (Lahn and Page 1999; Zhou 185 et al. 2014). While Amh is probably the candidate Z-linked male-determining gene given its 186 highly conserved function across vertebrates in testis development (Adolfi et al. 2019), we found 187 few expressed genes located within the S0 region of chrW except for one gene TIMM44. This 188 gene is essential for mitochondrial functions (Bonora et al. 2006), but has not been reported to 189 be participating in any sex determination process in other vertebrate species. Its chicken 190 ortholog has a gonad biased expression, and its Vac Z-linked copy has a biased expression in 191 the male brain (Supplementary Fig. 3), and the W-linked copy has a biased expression in the 192 ovary (Figure 2d). Whether this gene is indeed a female determining gene of Vac requires 193 future studies, particularly using the embryonic gonad samples when sex is determined. 194 195 7 10) o d .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 8 196 Figure 2 Evolutionary strata and candidate sex determining genes of Varanus acanthurus. a) 197 Vac sex chromosomes consist of at least three evolutionary strata. From top to bottom, we 198 show the male-to-female ratio of DNA read coverage; pairwise sequence identity between the 199 chrZ and chrW; and the heterozygosity level (SNP density) calculated by male reads along the 200 chrZ. b) Donut plot illustrates the degrees of degeneration of W-linked genes, divided by their 201 residing evolutionary strata (S1 and S0). The inner ring shows the genes with gametologs 202 present on both chrZ and chrW (‘with gametologs’), among which chrW gametologs are divided 203 into those with intact ORFs and robust gene expression (‘functional’, green), with disrupted 204 ORFs (grey), and with no expression (black). ‘without gametolog’ refers to genes whose chrW 205 gametologs have been deleted. c) Boxplots show the pairwise dN/dS values of 1:1 orthologs 206 between Phrynosoma platyrhinos (Ppl) and Varanus acanthurus (Vac), divided into those 207 located on the autosomes of both species (A-A), those on the autosome in Ppl, but in S0 208 stratum of chrZ in Vac (A-ZS0), and those on autosome in Ppl, but S1 of chrZ in Vac (A-ZS1). d) 209 The heatmap displays the z-scaled normalized expression levels (TPM) of S0 genes on the 210 chrZ, which were reported to be involved in sex-determination pathways in other vertebrates 211 (see Supplementary Table S3), along with all expressed chrW genes and their chicken 212 orthologs across different tissues. 213 214 Stepwise propagation of transposable elements in varanids 215 An enrichment of TEs at both chrW and the end of chr2 (termed ‘2q-end’ hereafter, region 216 defined by the sharp change of overall repeat content along chr2, Methods) (Figure 1c, 217 Supplementary Fig. 1) supports one of our hypotheses that the reported cytogenetic homology 218 is contributed by shared accumulation of TEs in both regions. We next ask when did such a TE 219 burst occur in both regions and what are the specific TEs that accumulated. To address this, we 220 first compared the overall TE content in the homologous regions of 2q-end of Vac, another 221 anguimorph lizard Elgaria multicarinata (Emu), and the Anguimorpha outgroup Ppl. An 222 enrichment of TEs, but to a significantly lower degree (P < 0.05, Wilcoxon test), with a much 223 shorter spanned region (23 Mb vs. 47 Mb of Vac), was found at the 2q-end of Ppl compared to 224 the other two species. This suggests that the 2q-end was already heterochromatic, and 225 probably has a low recombination rate in the ancestor of Anguimorpha dated 157 million years 226 ago (Figure 3a, Supplementary Fig. 4). Vac, Emu and Ppl together exhibit a gradient of TE 227 enrichment level and influenced region length at the 2q-end, supporting that this region 228 experienced stepwise propagation of TEs from the ancestor of Anguimorpha to that of varanids. 229 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 9 Particularly, we identified a long interspersed nuclear element (LINE) retrotransposon 230 family RTE-RovB that is enriched in both chrW and 2q-end compared to the rest of genome 231 (including chr4, whose end is also repetitive, Figure 1c, Supplementary Fig. 5) in both Vac 232 and another varanid V. salvator (Vsa), but not in the outgroup lizard species (Figure 3b). While 233 other retrotransposons, such as the long terminal repeat (LTR) families Gypsy and ERV1, show 234 different enrichment levels between 2q-end and chrW in the two varanids. An overall 14 fold 235 higher genome-wide percentage of RTE-RovB in varanids than their outgroup Anguimorpha 236 species suggests an ancestral burst of this TE family in the varanid ancestor, which then 237 preferentially inserted into the chrW and 2q-end with no or low recombination, consistent with 238 the "TE refugium hypothesis" (Peona et al. 2021). It is also likely that the RTE-RovB elements 239 inserted into either chromosomal region followed by transpositions onto another. 240 The TE burst may not be the only cause of the cytogenetic homology. And a more 241 functionally relevant consequence is that it can mediate interchromosomal segmental 242 duplications between the chrW and 2q-end. To investigate this, we aligned the genome 243 sequences of both chrZ and chrW against the chr2, after masking the repeat content. The 244 alignment of chrZ is supposed to inform us of duplications between chrW and 2q-end, followed 245 by loss of genes on the chrW. A total of 45 paralogs of Z-linked but not W-linked genes was 246 found on the short arm of chr2 (2p) but not in the 2q-end, 73% of which are shared with 247 chickens (Supplementary Fig. 6), suggesting ancient gene duplications predating the 248 divergence between birds and reptiles. We showed previously that the chicken homologous 249 chromosomes chrZ and chr28 of the Vac chr2p and chrZ (Figure 1d) are derived from the same 250 ancestral vertebrate chromosome approximated by that of amphioxus (Huang et al. 2023). Here 251 the duplicated genes shared between the Vac chr2p and chrZ are likely products of whole 252 genome duplication dated to the vertebrate ancestor, commonly referred to as ohnologs. Indeed, 253 among all duplicated gene pairs on the chr2p and chrZ, 67% of them have been previously 254 annotated as ohnologs (Singh et al. 2015). 255 By contrast, the non-repetitive homologous sequences between chr2 and chrW are only 256 found at the 2q-end, and are concentrated at the last 20Mb region (Figure 3c). The RTE-BovB 257 elements but not other TEs are also significantly (P < 0.05, Wilcoxon test) enriched within this 258 20Mb region. These results suggest that there could be one large segmental duplication (the 259 last 20Mb region) mediated by and inserted with the RTE-BovB elements, followed by other 260 individual duplications elsewhere in the 2q-end. We found 274 genes of chr2 with one or more 261 duplicated copies on the chrW, 90% of which are located within the 2q-end. If gene duplications 262 between chr2q and chrW were mediated by the varanid burst of RTE-BovBs, we expected to 263 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 10 find enrichment of this TE family nearby the duplicated genes. Indeed, among all subfamilies of 264 RTE-BovB, we identified md-5_family-305 and md-5_familiy-1095 that are much more enriched 265 nearby the paralogs between the 2q-end and chrW compared to other single copy genes as a 266 control. Other types of TEs are enriched to a similar level nearby the paralogs and the single 267 copy genes (Figure 3d). In addition, we performed Transpositions in Transpositions (TinT) 268 analyses of nested TEs, which assumes the younger TEs inserted more likely into the older 269 ones than the opposite direction of insertion. The result supported that these two RTE-BovB 270 subfamilies are indeed younger than other subfamilies (Figure 3e). Consistent with this, copies 271 of these two subfamilies located at the 2q-end and chrW have a significantly (P < 0.05, 272 Wilcoxon test) longer length, i.e., are more intact, but a significantly (P < 0.05, Wilcoxon test) 273 lower divergence level from their consensus sequences, than those elsewhere in the genome 274 (Figure 3f, Supplementary Fig. 7). Finally, when looking into the transcription of TEs across all 275 tissues, we found only the md-5_family-1095 of RTE-BovB have pronounced expression, 276 indicating its activity (Figure 3g). All these results showed that due to the no or low 277 recombination in chrW and 2q-end regions, certain active and evolutionarily young RTE-BovB 278 subfamilies experienced a copy number burst and preferentially accumulated at both regions in 279 the varanid ancestor, which probably in turn mediated segmental and gene duplications 280 between the two regions. 281 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 11 282 Figure 3 Gene duplications mediated by transposable elements between chrW and chr2q. a) 283 Chromosome alignments of the sex chromosomes and chr2 in Varanus acanthurus (Vac) 284 compared to another anguimorph lizard, Elgaria multicarinata (Emu), and the outgroup species 285 Phrynosoma platyrhinos (Ppl). The black lines within each homologous chromosome of different 286 species represent the density of overall TEs that shows a propagation at the 2q-end ( Figure 1c)287 Red barcodes below indicate the number (scaled to the red color) of duplicated genes per 100 288 kb genomic window between chr2 and chrW. b) Heatmap compares the normalized repeat 289 abundance between chrW, 2q-end, other chr2 regions, and other autosomes across 290 representative anguimorph lizards and one outgroup lizard Ppl. Only LINE/RTE-BovB is 291 enriched at the distal end of varanus lizards (Vac and Vsa) 2q-end and the chrW, compared to 292 other genomic regions. c) Diagrams of the zoomed-in 2q-end region, with the first bar plot 293 showing aligned unique region lengths with chrW. Below are line plots showing the lengths of 294 TEs in 100kb sliding windows. The region with grey background in the first bar plot corresponds 295 11 nt ). s .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 12 to areas with higher gene density in the barcode. d) Bubble plots show the enriched levels of 296 each TE subfamily nearby single-copy genes (dots) or duplicated genes (triangles), and those of 297 genes near each TE subfamily (Methods). e) Inferred relative ages of the subfamilies of RTE-298 BovB by TinT analysis, indicated by ovals and lines representing relative time points. f) Boxplots 299 show the distribution of Kimura distance, another measurement of TE relative age, of all copies 300 of rnd-5_family-1095 and rnd-5_family-305 subfamilies across different genomic regions (Macro 301 refers to all macro chromosomes except 2q-end region and micro refers to all micro 302 chromosomes except chrW). g) Heatmap shows the normalized expression levels of all RTE-303 BovB subfamilies across all tissues. 304 305 Gene traffic shaped the evolution of varanid sex chromosomes 306 It still remains unclear whether the varanid TE burst had mediated the genes duplicated from 307 2q-end to chrW or the reverse. This could respectively reflect a different evolutionary force of 308 either female-specific selection (Moghadam et al. 2012) that recruit genes elsewhere, or rescue 309 of degenerating genes, on the chrW (Hughes et al. 2015). To discern the two scenarios, we first 310 identified four multicopy gene families that have members located on both 2q-end and chrW in 311 Vac, and also have a homolog with well annotated function in mouse. The phylogenetic trees 312 constructed with these Vac gene sequences and their homologs in other non-varanid species, 313 which experienced a lower degree of TE propagations (Figure 3a), thus likely have fewer or no 314 duplications, are expected to reveal the direction of duplications into or out of chrW. These four 315 gene families consist of V2R, RNF39, HLA-DPB1, and EEF2, with a total of 442 gene copy 316 numbers mainly distributed on chr2 and the Z and W sex chromosomes (Figure 4a). All but 317 EEF2 have a significantly (Chi-squared test, P < 0.05) higher copy number in Vac than in three 318 other non-varanid lizards (Emu, Ppl, and Tiliqua scincoides or Tsc). Both RNF39 and HLA-319 DPB1 are associated with immune regulation(Liu et al. 2021; Wang et al. 2021), and V2R or the 320 vomeronasal chemosensory Type 2 receptor gene family encodes the pheromone receptors for 321 environmental perception that are associated with mate choice and predator avoidance(Silva 322 and Antunes 2017). 323 V2R, RNF39, HLA-DPB1 are mainly concentrated on the 2q-end as large tandem gene 324 clusters, but are randomly distributed along the chrW and chrZ, with V2R having additional large 325 numbers of copies on one end of chr4 (Figure 4b). Among the three families, the V2R family 326 has probably experienced the most pronounced expansion of copy numbers in the ancestor of 327 varanids. We annotated a total of 311 V2R copies throughout the genome of Vac, and both Vac 328 and Vsa hav e about 4-fold and more than 10-fold more gene copies on the 2q-end than the 329 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 13 outgroup species Emu, and even larger (> 10-fold) than in Ppl and Tsc (Figure 4c). While the 330 V2Rs on chr4 have much less copies in the varanids than in the other lizards. Our gene trees 331 support a different evolutionary history of V2Rs on different chromosomes: the V2Rs of chr4 332 form separate lineages from those of chr2 and chrW that are clustered with each other. Two 333 lineages are of particular interests that can inform the direction of duplications. One lineage 334 includes only the homologous V2Rs on the chrZ of Emu (whose chrW sequence is not 335 available), those on the homologous autosome of Tsc, both of which are then clustered with 336 copies from chrZ, chrW and chr2 of Vac. This led to an estimate of 17 V2R copies located on 337 the 2q-end of Vac as duplication products from chrZ or chrW. While another lineage clusters the 338 V2Rs of both chrW and chr2 of varanids with the chr2 copies of outgroup species, suggesting 339 the reverse direction of duplications of 19 W-linked copies from chr2 in the Vac and Vsc (Figure 340 4d). We thus found duplications of both directions, i.e., evidence of gene traffic between the sex 341 chromosomes and chr2 from V2Rs. For HLA-DPB1 and RNF39, gene duplications seem to be 342 exclusively from the 2q-end to chrW, generating only one gene copy onto the chrW. And EEF2 343 probably do not have interchromosomal duplications, according to their gene trees with other 344 lizards (Supplementary Fig. 8). 345 In particular, there are 11 V2R copies on the chrW, and they comprise the largest W-346 linked multigene family in Vac that is also female specific (Supplementary Fig. 9). As expected, 347 4 or 36% of these V2Rs have a disrupted ORF, due to either premature stop codon or 348 frameshift mutations. This percentage is much higher than that of V2Rs on the autosomes, 349 suggesting some of these W-linked copies are degenerating (Supplementary Fig. 10). Indeed, 350 the evolution rates (dN/dS ratios) of V2R copies without an intact ORF are significantly (P < 351 0.05, Wilcoxon test) faster than those with an intact ORF, on both chrW and autosomes 352 (S upplementary Fig. 11). Interestingly, all of the Vsa W-linked V2Rs (10 copies) have an intact 353 ORF (Supplementary Fig. 10). In addition, when we compare the evolutionary rates (pairwise 354 dN/dS ratios calculated with homologs of Emu) of V2Rs to other multicopy gene families whose 355 copy numbers are similarly abundant in the Vac genome, all the V2R copies, including those on 356 the chrW, show a significantly lower evolutionary rate than other multicopy gene families 357 (Figure 4e), indicating some selective constraints. All these results suggest that despite being 358 under the degenerating genetic environment of chrW, some V2R copies nevertheless have 359 retained functions after transposition and propagation. 360 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 14 361 Figure 4 Gene traffic between the varanid chrW and chr2. a) Gene duplications between Vac 362 Chr2 and sex chromosomes mainly belong to four gene families: V2R, RNF39, HLA-DPB1 and 363 EEF2, with the bubble size scaled to the copy number on each chromosome. b) Barcodes show 364 the genomic positions of V2R, RNF39, HLA-DPB1 and EEF2 family members in the Vac 365 genome, which are mainly concentrated at the chromosome ends of chr2 and chr4. c) Barplot 366 depicted the copy numbers of V2R genes in several lizard species, including Vac (V. 367 acanthurus), Vsa (V. salvator), Emu (E. multicarinata), Ppl (P. platyrhinos), and Tsc (Tiliqua 368 scincoides), with different color showing the copies from different chromosomes. d) A 369 maximum-likelihood tree constructed using the protein sequences of V2Rs from Vac, Vsa, Emu, 370 Ppl and Tsc (different colors represent different species in the inner ring), with the main 371 branches shown by different background colors as grey or white. W-linked copies (red color in 372 the outter ring) are mainly clustered with chr2 copies (orange color). e) Boxplots show the 373 pairwise dN/dS values between 1:1 orthologs in Emu and Vac, which are divided into single 374 copy genes, the top three multicopy (copy number > 3) gene families with most abundant copy 375 numbers throughout the genome, including HLA-DPB1, MR1 and RNF39 as a control, and 376 V2Rs from different chromosomes. 377 378 14 w , .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 15

379 Y or W chromosome regions with complete suppression of recombination are expected to 380 accumulate massive transposable elements and suffer long-term functional deterioration and 381 loss of gene functions(Charlesworth and Charlesworth 2000). Therefore, chrY/W is expected to 382 be an unfavored destination for gene transposition. In fact, independent Y-to-autosome 383 transpositions have been reported in mammals, probably as an escape from the degenerating 384 chrY to rescue the gene functions(Hughes et al. 2015). Contrary to this assumption, individual 385 cases of gene acquisition on the chrY/W have been reported in Drosophila (Koerich et al. 2008; 386 Carvalho et al. 2015; Tobler et al. 2017; Ricchio et al. 2021), in mammals (Page et al. 1984; 387 Murphy et al. 2006; Li et al. 2013; Janeč ka et al. 2018) and in birds (Xu et al. 2020), suggesting 388 the evolutionary course of these sex specific chromosomes can be more complex than the 389 canonical trajectory of degeneration. Our work here indicates one of such scenarios that is 390 mediated by a recent and local burst of TEs on an autosome, followed by directional 391 transpositions of a large number of genes onto the chrW. 392 Apart from the telomeric/centromeric or nucleolar organizer regions, previous works 393 have reported many cases of regional burst of certain TEs on one pair of autosomes (Oliveira 394 and Wright 1998), or one chromosome of an autosome pair (Bezerra et al. 2012; Moreira et al. 395 2013), with the latter case often confounding the identification of heteromorphic sex 396 chromosomes. For example, the long arm of one chr6 of platypus is enriched for satellite 397 repeats and ribosomal genes (Rens et al. 2004). Our recent work identified some isolated 398 populations of Vac and V. citrinus have an autosome showing cytogenetic homology to the 399 entire chrW (Dobry et al. 2025). While over one third of the entire chr3 of several tilapia species 400 consist of heterochromatin, which could predispose the region for becoming a sex determining 401 region (Tao et al. 2021). Other regional repeat expansion identified in humans frequently 402 underlie many diseases, including amyotrophic lateral sclerosis (ALS) and the fragile X 403 syndrome (Malik et al. 2021). Although the molecular and evolutionary causes of such local TE 404 bursts at autosomal regions are unclear, the sex specific chromosomes or other genomic 405 regions of low recombination, like the 2q-end identified in this work, due to their inefficiency of 406 resisting the TE transpositions, are more likely to be influenced by the burst than the rest 407 recombining regions. When the microchromosomes or other gene-rich regions undergo such TE 408 bursts, a byproduct of TE transposition onto the chrY/W is acquisition of a large number of 409 genes from the autosomes like what we observed in Vac (Figure 5). 410 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 16 411 Figure 5 Model of gene expansion mediated by transposable elements (TEs) on sex 412 chromosomes of Varanus lizards. This figure shows how consecutive expansions of TEs on 413 autosomes during speciation of anguimorph lizards had mediated expansion of genes on their 414 shared chrW. We inferred the first burst of TEs on the end of chr2 occurred in the Anguimorpha 415 ancestor, and the second in the Varanus ancestor. The second burst has likely caused 416 transpositions of TEs and their nearby genes between the chrW and chr2, which results in the 417 cytogenetic homology shared only between chrW and chr2 observed in Figure 1. 418 The genes transposed onto the chrY/W, are then expected to undergo functional 419 degeneration in the long-term, but at the same time, are subjected to immediate sex-specific 420 selection. For example, in both Drosophila and many mammals including human, cat and cattle, 421 some autosome-derived Y-linked single-copy gene (Carvalho et al. 2015) or multigene families 422 (Murphy et al. 2006; Yang et al. 2011; Chang et al. 2013) have acquired testis-specific 423 expression, in contrast to their female-biased or unbiased autosomal progenitors. Some 424 autosome derived gene fragments were also found to be amplified on the W chromosome of a 425 lizard Eremias velox (Lisachov et al. 2021) and Australian dragon lizard Pogona vitticeps (Ezaz 426 et al. 2013; Matsubara et al. 2019; Alam et al. 2020). An intensively characterized example is 427 the Y-linked DAZ gene family associated with human male fertility, which originated from the 428 single copy autosomal DAZL gene and expanded the copy number on the Y chromosome 429 16 a e, .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 17 during the primate evolution (Saxena et al. 1996; Yu et al. 2008). The example of V2R that we 430 reported in this work (Figure 5) represents the largest autosome-to-sex chromosome 431 transposition known to date. A previous study of the draft genome of another varanid lizard 432 Komodo dragon (Varanus komodoensis) reported expansion of V2R (Lind et al. 2019) and here 433 we concluded they are derived from duplications between chr2 and chrW, mediated by the TE 434 burst. The shared large number of intact W-linked V2Rs between varanids during the 28.2 435 million years of evolution, and the pattern of evolutionary rates suggest many of them have 436 retained functions (Figure 4), and possibly contribute to the sexually dimorphic perception of 437 environmental cues in Vac. However, further investigation is required to demonstrate their 438 detailed functionality, particularly tissue specific expression profile in those organs associated 439 with expression of V2Rs, including the vomeronasal organ. 440 441

442 Sample collection 443 Individuals used in the study were collected under authorization of the University of Canberra 444 Animal Ethics committee (Project ID: 20180306), and collection permits were obtained from the 445 Northern Territory (Permit number 63414) and Queensland (Permit number WA0010049) 446 governments for scientific collection. They were imported into the ACT under an import license 447 (LT201829). All methods reported in this study were conducted in accordance with the relevant 448 guidelines and regulations. The study is reported in accordance with ARRIVE guidelines. 449 Collection details have been described by Dobry et al, in previous studies (Dobry et al. 2023a; 450 Dobry et al. 2023b). Individuals were air freighted to the Canberra Airport, transported to the 451 University of Canberra, and housed in terrariums as described by Retes and Bennett (Daniel 452 and Frank 2001). Live tail tissue was used to culture cells for cytogenetic analysis as described 453 previously (Dobry et al. 2023a; Dobry et al. 2023b). Euthanasia was carried out using 454 intraperitoneal injection of sodium pentobarbitone at 60mg/kg sodium pentobarbital. Tissues 455 were immediately collected and snap frozen with liquid N2 and then processed for transcriptome 456 and DNA extraction. 457 458 Cytogenetic works 459 Animal collection, microdissection, preparation of sex chromosome specific probes and 460 validation of probes were described in our previous studies (Dobry et al. 2023b; Dobry et al. 461 2025). FISH mapping of sex chromosome probes are described in (Dobry et al. 2025). Briefly, 462 sex chromosome probe sequences composed of short, labelled polynucleotides ranging from 463 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 18 45-47 nucleotides in length were custom designed and manufactured by Arbor Biosciences 464 (myTags, Arbor Biosciences, Ann Arbor, MI, USA). These probes were designed from single 465 copy regions of the sex specific chromosome scaffolds identified in a previous study (Zhu et al. 466 2022) and validated in (Dobry et al. 2025). Probes were used at a concentration of 100 ng/μ L, 467 with a total of 200 ng per slide with 38 μ L hybridization buffer (BioCare Medical, Pacheco, CA). 468 Coverslips were placed on each slide and sealed with rubber cement to avoid dehydration 469 during denaturing at 68 °C for 5 min and then incubated at 37 °C for 48 h. The coverslips were 470 then removed and the slides were washed with 0.4x SSC (3 M NaCl, 0.3 M sodium citrate, pH 7, 471 and 0.3% (v/v) IGEPAL (Sigma-Aldrich)) at 60 °C for 2 min. A second wash was then performed 472 with 2x SSC for 1 min. The slides were then dehydrated with a series 1 min ethanol washes 473 consisting of 70%, 90%, and 100% (v/v) respectively and allowed to completely dry before 474 staining with Vectashield antifade mounting medium with DAPI (Vector Laboratories). Slides 475 were viewed and photographed with a Leica Microsystems Thunder Imaging system, and then 476 karyotype images were constructed using Adobe Photoshop 2023. 477 478 Genome assembly and annotation 479 We sequenced and assembled the reference genome of a female V. acanthurus. Liver tissues 480 were collected for DNA extraction and library preparation for Illumina, PacBio SMRT, and Hi-C 481 sequencing (See Supplementary Table S1). Original assembly was performed with Flye (v2.7-482 b1585) (Kolmogorov et al. 2019) using PacBio reads and polished with Illumina DNA reads 483 using Pilon (v1.22) (Walker et al. 2014). The assembly was subsequently scaffolded with 484 Ragtag (v2.1.0) (Alonge et al. 2022) using BGI stLFR reads, and completed with the 3D-DNA 485 pipeline (v180922) (Dudchenko et al. 2017) based on Hi-C data, resulting in an assembly with 486 over 97% of sequences assigned to chromosomes, and BlobToolKit (v4.3.10) (Challis et al. 487 2020) were later used for genome quality control. For genome annotation, a de novo repeat 488 sequence database was first constructed using RepeatModeler (v2.0.3) (Flynn et al. 2020), and 489 combined with the squamate repeat consensus library from Repbase (v20181026) (Bao et al. 490 2015) to create a comprehensive repeat library. Repeat sequences in the genome were then 491 identified and classified with RepeatMasker(Smit et al.) (v4.07). Later we used the Funannotate 492 pipeline (v1.8.5) (Palmer and Stajich) for the gene annotation, obtaining transcriptome data from 493 brain and gonad tissues, each with two replicates per sex for comprehensive transcript 494 assembly. To construct a non-redundant reference protein library, we included proteins from 495 Komodo dragon (Varanus komodoensis, GCF_004798865.1, (Lind et al. 2019)), green anole 496 (Anolis carolinensis, GCF_000090745.2, (Alföldi et al. 2011)), Indian cobra (Naja naja, 497 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 19 GCA_009733165.1, (Suryamohan et al. 2020)), chicken (Gallus gallus, GCF_016699485.2, 498 (Rhie et al. 2021)), mouse (Mus musculus, GCF_000001635.27, (Church et al. 2011)), and 499 human (Homo sapiens, GCF_000001405.40, (Schneider et al. 2017)) using cd-hit(v4.7) (Fu et 500 al. 2012). Among all annotated genes, those with transcripts containing more than two 501 premature stop codons or frameshift mutations were classified as having disrupted open 502 reading frames (ORFs), and genes with more than two copies were classified as multicopy 503 genes. 504 505 Sex chromosome identification 506 To identify the sex chromosome sequences of V. acanthurus, Illumina DNA reads from both 507 sexes were aligned to the masked genome using Bowtie2 (v2.2.9) (Langmead et al. 2019). 508 Read coverage for each sex was calculated in 100-kb non-overlapping windows using SAMtools 509 (v1.6) (Danecek et al. 2021) and normalized to the median depth per base pair across the entire 510 genome. This normalization facilitated direct comparison of coverage between sexes. The 511 "covered site" was defined as a base pair with a read coverage of at least 3. The Z chromosome 512 was expected to display an autosome-like male-to-female (M/F) coverage ratio but with half the 513 overall read coverage due to hemizygosity in females. Conversely, the highly differentiated W 514 chromosome was identified by scaffolds exhibiting both an M/F coverage ratio and a covered 515 site ratio below 0.2, alongside a scaffold size exceeding 5 kb. These W-linked scaffolds were 516 subsequently assembled into a contiguous W chromosome using Ragtag (v2.1.0), with the Z 517 chromosome serving as a reference (Supplementary Table S4). For V. salvator, where 518 Illumina DNA reads (SRR16080542) were available only from female tissues, autosomes and 519 sex-linked sequences were distinguished based on sequencing depth. Autosomal sequences 520 displayed approximately two-fold the depth of sex-linked sequences. Homology-based 521 comparisons with the Z and W chromosomes of V. acanthurus were then employed to classify 522 W-linked and Z-linked sequences in V. salvator. 523 524 Evolutionary strata 525 To identify the evolutionary strata of the sex chromosomes in Varanus acanthurus, we first 526 identified regions with similar sequencing depth between sexes, as the pseudoautosomal region 527 (PAR). The sequence similarity between the Z and W chromosomes was then assessed by 528 aligning the masked W chromosome sequence to the Z chromosome using Lastz (v1.02.00) 529 (Harris 2008) with parameters “--hspthresh=2200 --inner=2000 --ydrop=3400 --530 gappedthresh=10000”. Syntenic fragments were merged into longer blocks through the 531 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 20 alignment chains and nets by the UCSC pipeline. Sequence similarity between the sex 532 chromosomes was measured in 100-kb sliding windows along the Z chromosome. To evaluate 533 male heterozygosity along the Z chromosome, SNP calling was performed using the 534 HaplotypeCaller module from the Genome Analysis Toolkit (GATK, v3.8) (Van der Auwera and 535 O'Connor) with male Illumina DNA reads. SNPs were filtered under the following criteria: "QD < 536 2.0, MQ 60.0, SOR > 3.0, MQRankSum < −12.5, and ReadPosRankSum < −8.0." 537 SNP counts for each chromosome were calculated in 100-kb windows. Both sequence similarity 538 and male heterozygosity values were used to change-point analysis using the R package cpm 539 (v2.3) (Ross 2015). The resulting strata were statistically evaluated using the Wilcoxon-test to 540 confirm significant differences between strata. 541 542 Gene expressions 543 To quantify the gene expressions, RNA-seq reads from brain tissues of Varanus acanthurus of 544 both sexes, as well as testis and ovary, were mapped to the genome using Hisat2 (v2.1.0) (Kim 545 et al. 2019), retaining only uniquely mapped reads (flagged as NH:i:1 and without the ZS:i flag). 546 FeatureCounts (v1.6.2) (Liao et al. 2014) were used to count these reads, and gene expression 547 levels were normalized and quantified as TPM (Transcripts Per Million). A similar pipeline was 548 employed to quantify gene expression levels in chicken, utilizing RNA-seq reads from liver and 549 muscle tissues of both sexes, as well as testis and ovary (see Supplementary Table S2). For 550 all figures, the expression levels of each gene in each tissue were represented by the median 551 values across replicates. Genes with TPM < 1 across all tissues were classified as silenced. 552 553 Inferring synteny between genomes 554 To identify synteny blocks between species in this study, reciprocal best-hit (RBH) blast was 555 performed using BLASTp (v2.6.0) to detect one-to-one orthologs based on the protein 556 sequences of all species. The search employed a stringent e-value threshold of 1e-10, and 557

were filtered to retain alignments covering more than 50% of the reference protein and 558 exhibiting at least 50% sequence identity. In addition to Varanus acanthurus, the species 559 included in this comparative study were another monitor lizard (Varanus salvator, NCBI 560 accession JAIXND000000000.1, (Chetruengchai et al. 2022)), a lizard from Anguimorpha 561 (Elgaria multicarinata, GCF_023053635.1), a lizard from Iguania (Phrynosoma platyrhinos, 562 GCA_020142125.1, (Koochekian et al. 2022)), a skink (Tiliqua scincoides, GCF_035046505.1, 563 (Rhie et al. 2021)), and chicken (Gallus gallus, GCF_016699485.2, (Rhie et al. 2021)). 564 565 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 21 dN/dS calculation 566 Pairwise dN/dS values between orthologs were calculated as follows: one-to-one orthologs 567 were first identified using the RBH pipeline described earlier. Coding sequences for each 568 ortholog pair were aligned using Prank (v150803) (Löytynoja 2014) with default parameters. 569 dN/dS values were then calculated for each pairwise alignment using the codeml module of 570 PAML (v4.8) (Yang 2007)under the free-ratio model. The same pipeline, with identical 571 parameters and modes, was applied to analyze orthologs between Varanus acanthurus and 572 Elgaria multicarinata, as well as Varanus acanthurus and Phrynosoma platyrhinos. 573 574 Repeat analysis 575 To assess the relative abundance of retroposons at the distal ends of chromosomes 2 and the 576 whole W in varanids compared to homologous genomic regions in the other species, we 577 calculated the length of each repeat family in 100kb windows across these chromosomes and 578 normalized by the chromosome-wide mean. To investigate the specific repeat subfamilies 579 potentially driving gene duplication events, repeats located within ±10kb of a gene were 580 categorized as nearby repeats. To quantify enrichment, we calculated the proportion of genes 581 with nearby repeats from each subfamily, along with the proportion of copy number of these 582 nearby repeats relative to total loci on the chromosomes. To assess the relative activity and age 583 of retroposons, we first applied TinT (Transpositions in Transpositions) algorism to examine the 584 frequencies of nested transpositions. We also used RepeatMasker outputs to calculate Kimura 585 distances for each retroposon subfamily to measure the divergence level of each element. For 586 further phylogenetic inference, consensus sequences of the RTE/BovB subfamilies from two 587 varanids and the desert horned lizard (Ppl) were analyzed using IQ-TREE. 588 RNA-seq data from brain tissues of both sexes, male testis, and female ovary were mapped to 589 the genome using the same pipeline as employed for gene expression analysis. The resulting 590 BAM files were subsequently processed with TEcount (v1.0.1) (Marasca et al. 2022) for 591 quantification of counts by class, family and subfamily levels. Repeat expression levels were 592 later normalized as CPM (Counts Per Million). For all figures, the expression levels in each 593 tissue were represented by the median values across replicates. 594 595 Paralog analysis 596 To first identify paralogs in the genome of each species, the longest transcript of each gene in 597 Vac, Vsa, Emu, Ppl and Tsc were fed to OrthoFinder (v2.5.4) (Emms and Kelly 2019) with 598 default parameters, to get orthogroups that pinpoint inter- and intra-species duplications. For 599 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 22 the identification of ohnologs, the defined duplications were queried against the database at 600 http://ohnologs.curie.fr/, using the chicken data as a reference. Later, protein sequences of all 601 paralogs within each orthogroup were first aligned using MAFFT (v7.505) (Katoh and Standley 602 2013) with default parameters, and poorly aligned regions were trimmed with TrimAl (v1.5.rev0) 603 (Capella-Gutiérrez et al. 2009) with -gt 0.1 option. Maximum likelihood trees were then inferred 604 using IQ-TREE (v1.6.11) (Minh et al. 2020). The direction of the duplication events in varanids 605 were inferred based on both the genomic positions of their homologs in the outgroup species 606 Emu, Ppl and Tsc, and the structure of phylogeny tree from each gene family. 607 608 Data Access 609 The genome assembly is available on NCBI under accession numberPRJNA855548. The 610 annotation file is available on Github: https://github.com/zjuzexian/Varanus-sex-chromosomes/. 611 Sequenced reads generated during this study are available on NCBI PRJNA1201623. The 612 stLFR reads, produced in our previous study, can be accessed under accession ID 613 PRJNA737594. Varanus acanthurus draft genome assembly that are used to develop sex 614 chromosome specific probes are available via genome warehouse PRJCA005583 and sex 615 chromosome probe sequences are available in Dobry et al (Dobry et al. 2025). Note the codes 616 used to generate the figures can be found on GitHub: https://github.com/zjuzexian/Varanus-sex-617 chromosomes/ 618 619

620 Qi Zhou is supported by the National Key Research and Development Program of China 621 (2023YFA1800500), National Natural Science Foundation of China (32170415). Jason Dobry is 622 partially supported by the Australian Government Research Training Program (RTP) stipend 623 scholarship. Tariq Ezaz and Erik Wapstra were partially supported by the Australian Research 624 Council Discovery Project grant (ARC DP200101406). 625 626

627 Adolfi MC, Nakajima RT, Nóbrega RH, Schartl M. 2019. Intersex, Hermaphroditism, and 628 Gonadal Plasticity in Vertebrates: Evolution of the Müllerian Duct and Amh/Amhr2 629 Signaling. Annu Rev Anim Biosci 7: 149-172. 630 Alam SMI, Altmanová M, Prasongmaneerut T, Georges A, Sarre SD, Nielsen SV, Gamble T, 631 Srikulnath K, Rovatsos M, Kratochvíl L et al. 2020. Cross-Species BAC Mapping 632 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 23 Highlights Conservation of Chromosome Synteny across Dragon Lizards (Squamata: 633 Agamidae). Genes 11: 698. 634 Alföldi J, Di Palma F, Grabherr M, Williams C, Kong L, Mauceli E, Russell P, Lowe CB, Glor RE, 635 Jaffe JD et al. 2011. The genome of the green anole lizard and a comparative analysis 636 with birds and mammals. Nature 477: 587-591. 637 Alonge M, Lebeigle L, Kirsche M, Jenike K, Ou S, Aganezov S, Wang X, Lippman ZB, Schatz 638 MC, Soyk S. 2022. Automated assembly scaffolding using RagTag elevates a new 639 tomato system for high-throughput genome editing. Genome Biology 23: 1-19. 640 Bachtrog D, Toda NRT, Lockton S. 2010. Dosage Compensation and Demasculinization of X 641 Chromosomes in Drosophila. Current Biology 20: 1476-1481. 642 Bao W, Kojima KK, Kohany O. 2015. Repbase Update, a database of repetitive elements in 643 eukaryotic genomes. Mobile DNA 6: 1-6. 644 Betrán E, Thornton K, Long M. 2002. Retroposed new genes out of the X in Drosophila. 645 Genome Res 12: 1854-1859. 646 Bezerra AMR, Pagnozzi JM, Carmignotto AP, Yonenaga-Yassuda Y, Rodrigues FHG. 2012. A 647 new karyotype for the spiny rat Clyomys laticeps (Thomas, 1909) (Rodentia, Echimyidae) 648 from Central Brazil. Comp Cytogenet 6: 153-161. 649 Bonora E, Evangelisti C, Bonichon F, Tallini G, Romeo G. 2006. Novel germline variants 650 identified in the inner mitochondrial membrane transporter TIMM44 and their role in 651 predisposition to oncocytic thyroid carcinomas. Br J Cancer 95: 1529-1536. 652 Bull JJ. 1980. Sex Determination in Reptiles. The Quarterly Review of Biology 55: 3-21. 653 Bull JJ, Charnov EL. 1985. ON IRREVERSIBLE EVOLUTION. Evolution 39: 1149-1155. 654 Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. 2009. trimAl: a tool for automated 655 alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25: 1972. 656 Carvalho AB, Vicoso B, Russo CAM, Swenor B, Clark AG. 2015. Birth of a new gene on the Y 657 chromosome of Drosophila melanogaster. Proc Natl Acad Sci U S A 112: 12450-12455. 658 Challis R, Richards E, Rajan J, Cochrane G, Blaxter M. 2020. BlobToolKit - Interactive Quality 659 As sessment of Genome Assemblies. G3 (Bethesda, Md) 10: 1361-1374. 660 Chang T-C, Yang Y, Retzel EF, Liu W-S. 2013. Male-specific region of the bovine Y 661 chromosome is gene rich with a high transcriptomic activity in testis development. Proc 662 Natl Acad Sci U S A 110: 12373-12378. 663 Charlesworth B, Charlesworth D. 2000. The degeneration of Y chromosomes. Philos Trans R 664 Soc Lond B Biol Sci 355: 1563-1572. 665 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 24 Charlesworth B, Coyne JA, Barton NH. 1987. The Relative Rates of Evolution of Sex 666 Chromosomes and Autosomes. The American Naturalist 130: 113-146. 667 Chetruengchai W, Singchat W, Srichomthong C, Assawapitaksakul A, Srikulnath K, Ahmad SF, 668 Phokaew C, Shotelersuk V. 2022. Genome of Varanus salvator macromaculatus (Asian 669 water monitor) reveals adaptations in the blood coagulation and innate immune system. 670 Front Ecol Evol 10: 850817. 671 Church DM, Schneider VA, Graves T, Auger K, Cunningham F, Bouk N, Chen HC, Agarwala R, 672 McLaren WM, Ritchie GR et al. 2011. Modernizing reference genome assemblies. PLoS 673 biology 9: e1001091. 674 Cortez D, Marin R, Toledo-Flores D, Froidevaux L, Liechti A, Waters PD, Grützner F, 675 Kaessmann H. 2014. Origins and functional evolution of Y chromosomes across 676 mammals. Nature 508: 488-493. 677 Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, 678 McCarthy SA, Davies RM et al. 2021. Twelve years of SAMtools and BCFtools. 679 GigaScience 10: giab008. 680 Daniel B, Frank R. 2001. Multiple generations, multiple clutches, and early maturity in four 681 species of monitor lizards (Varanidae) bred in captivity. Herpetological Review 32: 244-682 255. 683 Davalos-Dehullu E, Baty SM, Fisher RN, Scott PA, Dolby GA, Munguia-Vega A, Cortez D. 2023. 684 Chromosome-Level Genome Assembly of the Blacktail Brush Lizard, Urosaurus 685 nigricaudus, Reveals Dosage Compensation in an Endemic Lizard. Genome Biol Evol 686 15: evad210. 687 Dobry J, Wapstra E, Stringer EJ, Gruber B, Deakin JE, Ezaz T. 2023a. Widespread 688 chromosomal rearrangements preceded genetic divergence in a monitor lizard, Varanus 689 acanthurus (Varanidae). Chromosome research 31: 9. 690 Dobry J, Zhu Z, Zhou Q, Wapstra E, Deakin JE, Ezaz T. 2023b. Fixed Allele Differences 691 Associated With the Centromere Reveal Chromosome Morphology and Rearrangements 692 in a Reptile (Varanus acanthurus BOULENGER). Molecular biology and evolution 40: 693 msad124. 694 Dobry J, Zhu Z, Zhou Q, Wapstra E, Deakin JE, Ezaz T. 2025. The role of unbalanced 695 segmental duplication in sex chromosome evolution in Australian ridge-tailed goannas. 696 doi:10.21203/rs.3.rs-5547174/v1. 697 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 25 Dudchenko O, Batra SS, Omer AD, Nyquist SK, Hoeger M, Durand NC, Shamim MS, Machol I, 698 Lander ES, Aiden AP et al. 2017. De novo assembly of the Aedes aegypti genome using 699 Hi-C yields chromosome-length scaffolds. Science 356: 92-95. 700 Emerson JJ, Kaessmann H, Betrán E, Long M. 2004. Extensive gene traffic on the mammalian 701 X chromosome. Science 303: 537-540. 702 Emms DM, Kelly S. 2019. OrthoFinder: phylogenetic orthology inference for comparative 703 genomics. Genome Biology 20: 1-14. 704 Ezaz T, Azad B, O'Meally D, Young MJ, Matsubara K, Edwards MJ, Zhang X, Holleley CE, 705 Deakin JE, Marshall Graves JA et al. 2013. Sequence and gene content of a large 706 fragment of a lizard sex chromosome and evaluation of candidate sex differentiating 707 gene R-spondin 1. BMC Genomics 14: 1-13. 708 Ezaz T, Sarre SD, O'Meally D, Graves JAM, Georges A. 2009. Sex chromosome evolution in 709 lizards: independent origins and rapid transitions. Cytogenet Genome Res 127: 249-260. 710 Flynn JM, Hubley R, Goubert C, Rosen J, Clark AG, Feschotte C, Smit AF. 2020. 711 RepeatModeler2 for automated genomic discovery of transposable element families. 712 Proceedings of the National Academy of Sciences 117: 9451-9457. 713 Fu L, Niu B, Zhu Z, Wu S, Li W. 2012. CD-HIT: accelerated for clustering the next-generation 714 sequencing data. Bioinformatics 28: 3150. 715 Gamble T, Coryell J, Ezaz T, Lynch J, Scantlebury DP, Zarkower D. 2015. Restriction Site-716 Associated DNA Sequencing (RAD-seq) Reveals an Extraordinary Number of 717 Transitions among Gecko Sex-Determining Systems. Mol Biol Evol 32: 1296-1309. 718 Geneva AJ, Park S, Bock DG, de Mello PLH, Sarigol F, Tollis M, Donihue CM, Reynolds RG, 719 Feiner N, Rasys AM et al. 2022. Chromosome-scale genome assembly of the brown 720 anole (Anolis sagrei), an emerging model species. Commun Biol 5: 1126. 721 Harris RS. 2008. Improved pairwise alignment of genomic DNA. The Pennsylvania State 722 University. 723 Hill PL, Burridge CP, Ezaz T, Wapstra E. 2018. Conservation of Sex-Linked Markers among 724 Conspecific Populations of a Viviparous Skink, Niveoscincus ocellatus, Exhibiting 725 Genetic and Temperature-Dependent Sex Determination. Genome Biol Evol 10: 1079-726 1087. 727 Huang Z , Xu L, Cai C, Zhou Y, Liu J, Xu Z, Zhu Z, Kang W, Cen W, Pei S et al. 2023. Three 728 amphioxus reference genomes reveal gene and chromosome evolution of chordates. 729 Proc Natl Acad Sci U S A 120: e2201504120. 730 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 26 Hughes JF, Skaletsky H, Koutseva N, Pyntikova T, Page DC. 2015. Sex chromosome-to-731 autosome transposition events counter Y-chromosome gene loss in mammals. Genome 732 Biol 16: 104. 733 Iannucci A, Altmanová M, Ciofi C, Ferguson-Smith M, Milan M, Pereira JC, Pether J, Rehák I, 734 Rovatsos M, Stanyon R et al. 2019. Conserved sex chromosomes and karyotype 735 evolution in monitor lizards (Varanidae). Heredity 123: 215-227. 736 Janeč ka JE, Davis BW, Ghosh S, Paria N, Das PJ, Orlando L, Schubert M, Nielsen MK, Stout 737 TAE, Brashear W et al. 2018. Horse Y chromosome assembly displays unique 738 evolutionary features and putative stallion fertility genes. Nat Commun 9: 2945. 739 Johnson PM, Altmanová M, Rovatsos M, Velenský P, Vodič ka R, Rehák I, Kratochvíl L. 2016. 740 First Description of the Karyotype and Sex Chromosomes in the Komodo Dragon 741 (Varanus komodoensis). Cytogenetic and genome research 148: 284-291. 742 Katoh K, Standley DM. 2013. MAFFT Multiple Sequence Alignment Software Version 7: 743 Improvements in Performance and Usability. Molecular Biology and Evolution 30: 772. 744 Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. 2019. Graph-based genome alignment and 745 genotyping with HISAT2 and HISAT-genotype. Nature biotechnology 37: 907-915. 746 Kinga M, King D. 1975. Chromosomal Evolution in the Lizard Genus Varanus (Reptilia). Aust Jnl 747 Of Bio Sci 28: 89-108. 748 Koerich LB, Wang X, Clark AG, Carvalho AB. 2008. Low conservation of gene content in the 749 Drosophila Y chromosome. Nature 456: 949-951. 750 Kolmogorov M, Yuan J, Lin Y, Pevzner PA. 2019. Assembly of long, error-prone reads using 751 repeat graphs. Nature Biotechnology 37: 540-546. 752 Koochekian N, Ascanio A, Farleigh K, Card DC, Schield DR, Castoe TA, Jezkova T. 2022. A 753 chromosome-level genome assembly and annotation of the desert horned lizard, 754 Phrynosoma platyrhinos, provides insight into chromosomal rearrangements among 755 reptiles. Gigascience 11: giab098. 756 Krasovec M, Chester M, Ridout K, Filatov DA. 2018. The Mutation Rate and the Age of the Sex 757 Chromosomes in Silene latifolia. Curr Biol 28: 1832-1838.e1834. 758 Lahn BT, Page DC. 1999. Four evolutionary strata on the human X chromosome. Science 286: 759 964-967. 760 Langmead B, Wilks C, Antonescu V, Charles R. 2019. Scaling read aligners to hundreds of 761 threads on general-purpose processors. Bioinformatics (Oxford, England) 35: 421-432. 762 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 27 Leitão HG, Diedericks G, Broeckhoven C, Baeckens S, Svardal H. 2023. Chromosome-Level 763 Genome Assembly of the Cape Cliff Lizard (Hemicordylus capensis). Genome Biol Evol 764 15: evad001. 765 Li G, Davis BW, Raudsepp T, Pearks Wilkerson AJ, Mason VC, Ferguson-Smith M, O'Brien PC, 766 Waters PD, Murphy WJ. 2013. Comparative analysis of mammalian Y chromosomes 767 illuminates ancestral structure and lineage-specific evolution. Genome Res 23: 1486-768 1495. 769 Liao Y, Smyth GK, Shi W. 2014. featureCounts: an efficient general purpose program for 770 assigning sequence reads to genomic features. Bioinformatics (Oxford, England) 30: 771 923-930. 772 Lind AL, Lai YYY, Mostovoy Y, Holloway AK, Iannucci A, Mak ACY, Fondi M, Orlandini V, 773 Eckalbar WL, Milan M et al. 2019. Genome of the Komodo dragon reveals adaptations in 774 the cardiovascular and chemosensory systems of monitor lizards. Nat Ecol Evol 3: 1241-775 1252. 776 Lipinska AP, Toda NRT, Heesch S, Peters AF, Cock JM, Coelho SM. 2017. Multiple gene 777 movements into and out of haploid sex chromosomes. Genome Biol 18: 104. 778 Lisachov A, Andreyushkova D, Davletshina G, Prokopov D, Romanenko S, Galkina S, 779 Saifitdinova A, Simonov E, Borodin P, Trifonov V. 2021. Amplified Fragments of an 780 Autosome-Borne Gene Constitute a Significant Component of the W Sex Chromosome 781 of (Reptilia, Lacertidae). Genes (Basel) 12: 779. 782 Liu B, Shao Y, Fu R. 2021. Current research status of HLA in immune‐ related diseases. 783 Immunity, Inflammation and Disease 9: 340-350. 784 Löytynoja A. 2014. Phylogeny-aware alignment with PRANK. Multiple Sequence Alignment 785

doi:10.1007/978-1-62703-646-7_10: 155-170. 786 Malik I, Kelley CP, Wang ET, Todd PK. 2021. Molecular mechanisms underlying nucleotide 787 repeat expansion disorders. Nat Rev Mol Cell Biol 22: 589-607. 788 Mank JE, Nam K, Ellegren H. 2010. Faster-Z evolution is predominantly due to genetic drift. Mol 789 Biol Evol 27: 661-670. 790 Marasca F, Sinha S, Vadalà R, Polimeni B, Ranzani V, Paraboschi EM, Burattin FV, Ghilotti M, 791 Crosti M, Negri ML et al. 2022. LINE1 are spliced in non-canonical transcript variants to 792 regulate T cell quiescence and exhaustion. Nature Genetics 54: 180-193. 793 Matsubara K, O'Meally D, Sarre SD, Georges A, Srikulnath K, Ezaz T. 2019. ZW Sex 794 Chromosomes in Australian Dragon Lizards (Agamidae) Originated from a Combination 795 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 28 of Duplication and Translocation in the Nucleolar Organising Region. Genes (Basel) 10: 796 861. 797 Mezzasalma M, Guarino FM, Odierna G. 2021. Lizards as Model Organisms of Sex 798 Chromosome Evolution: What We Really Know from a Systematic Distribution of 799 Available Data? Genes 12: 1341. 800 Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams, von Haeseler A, Lanfear R. 801 2020. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the 802 Genomic Era. Molecular biology and evolution 37: 1530-1534. 803 Moghadam HK, Pointer MA, Wright AE, Berlin S, Mank JE. 2012. W chromosome expression 804 responds to female-specific selection. Proc Natl Acad Sci U S A 109: 8207-8211. 805 Moreira CN, Di-Nizo CB, de Jesus Silva MJ, Yonenaga-Yassuda Y, Ventura K. 2013. A 806 remarkable autosomal heteromorphism in Pseudoryzomys simplex 2n = 56; FN = 54-55 807 (Rodentia, Sigmodontinae). Genet Mol Biol 36: 201-206. 808 Murphy WJ, Pearks Wilkerson AJ, Raudsepp T, Agarwala R, Schäffer AA, Stanyon R, 809 Chowdhary BP. 2006. Novel gene acquisition on carnivore Y chromosomes. PLoS 810 Genet 2: e43. 811 Oliveira C, Wright JM. 1998. Molecular cytogenetic analysis of heterochromatin in the 812 chromosomes of tilapia, Oreochromis niloticus (Teleostei: Cichlidae). Chromosome Res 813 6: 205-211. 814 Page DC, Harper ME, Love J, Botstein D. 1984. Occurrence of a transposition from the X-815 chromosome long arm to the Y-chromosome short arm during human evolution. Nature 816 311: 119-123. 817 Palmer J, Stajich J. 2020. Funannotate v1.8.1: Eukaryotic genome annotation (v1.8.1). Zenodo 818 doi:https://doi.org/10.5281/zenodo.4054262. 819 Peona V, Palacios-Gimenez OM, Blommaert J, Liu J, Haryoko T, Jønsson KA, Irestedt M, Zhou 820 Q, Jern P, Suh A. 2021. The avian W chromosome is a refugium for endogenous 821 retroviruses with likely effects on female-biased mutational load and genetic 822 incompatibilities. Philos Trans R Soc Lond B Biol Sci 376: 20200186. 823 Pinto BJ, Nielsen SV, Sullivan KA, Behere A, Keating SE, van Schingen-Khan M, Nguyen TQ, 824 Ziegler T, Pramuk J, Wilson MA et al. 2024. It's a trap?! Escape from an ancient, 825 ancestral sex chromosome system and implication of Foxl2 as the putative primary sex-826 determining gene in a lizard (Anguimorpha; Shinisauridae). Evolution 78: 355-363. 827 Pokorná M, Kratochvíl L. 2009. Phylogeny of sex-determining mechanisms in squamate reptiles: 828 are sex chromosomes an evolutionary trap? Zool J Linn Soc 156: 168-183. 829 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 29 Rens W, Grützner F, O'Brien PCM, Fairclough H, Graves JAM, Ferguson-Smith MA. 2004. 830 Resolution and evolution of the duck-billed platypus karyotype with an 831 X1Y1X2Y2X3Y3X4Y4X5Y5 male sex chromosome constitution. Proc Natl Acad Sci U S 832 A 101: 16257-16261. 833 Rhie A McCarthy SA Fedrigo O Damas J Formenti G Koren S Uliano-Silva M Chow W 834 Fungtammasan A Kim J et al. 2021. Towards complete and error-free genome 835 assemblies of all vertebrate species. Nature 592: 737-746. 836 Ricchio J, Uno F, Carvalho AB. 2021. New Genes in the Y Chromosome: Lessons from D. 837 willistoni. Genes (Basel) 12: 1815. 838 Ross GJ. 2015. Parametric and Nonparametric Sequential Change Detection in R: The cpm 839 Package. J Stat Soft 66: 1-20. 840 Rovatsos M, Rehák I, Velenský P, Kratochvíl L. 2019. Shared Ancient Sex Chromosomes in 841 Varanids, Beaded Lizards, and Alligator Lizards. Mol Biol Evol 36: 1113-1120. 842 Sackton TB, Corbett-Detig RB, Nagaraju J, Vaishna L, Arunkumar KP, Hartl DL. 2014. Positive 843 selection drives faster-Z evolution in silkmoths. Evolution 68: 2331-2342. 844 Saxena R, Brown LG, Hawkins T, Alagappan RK, Skaletsky H, Reeve MP, Reijo R, Rozen S, 845 Dinulos MB, Disteche CM et al. 1996. The DAZ gene cluster on the human Y 846 chromosome arose from an autosomal gene that was transposed, repeatedly amplified 847 and pruned. Nat Genet 14: 292-299. 848 Schneider VA, Graves-Lindsay T, Howe K, Bouk N, Chen HC, Kitts PA, Murphy TD, Pruitt KD, 849 Thibaud-Nissen F, Albracht D et al. 2017. Evaluation of GRCh38 and de novo haploid 850 genome assemblies demonstrates the enduring quality of the reference assembly. 851 Genome research 27: 849-864. 852 Silva L, Antunes A. 2017. Vomeronasal Receptors in Vertebrates and the Evolution of 853 Pheromone Detection. Annu Rev Anim Biosci 5: 353-370. 854 Singh PP, Arora J, Isambert H. 2015. Identification of Ohnolog Genes Originating from Whole 855 Genome Duplication in Early Vertebrates, Based on Synteny Comparison across 856 Multiple Genomes. PLoS Comput Biol 11: e1004394. 857 Smit A, Hubley R, Green P. 2015. RepeatMasker Open-4.0. 2013–2015. Seattle, USA. 858 Sturgill D, Zhang Y, Parisi M, Oliver B. 2007. Demasculinization of X chromosomes in the 859 Dr osophila genus. Nature 450: 238-241. 860 Suryamohan K, Krishnankutty SP, Guillory J, Jevit M, Schröder MS, Wu M, Kuriakose B, 861 Mathew OK, Perumal RC, Koludarov I et al. 2020. The Indian cobra reference genome 862 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 30 and transcriptome enables comprehensive identification of venom toxins. Nature 863 Genetics 52: 106-117. 864 Tao W, Xu L, Zhao L, Zhu Z, Wu X, Min Q, Wang D, Zhou Q. 2021. High-quality chromosome-865 level genomes of two tilapia species reveal their evolution of repeat sequences and sex 866 chromosomes. Mol Ecol Resour 21: 543-560. 867 Tobler R, Nolte V, Schlötterer C. 2017. High rate of translocation-based gene birth on the Y 868 chromosome. Proc Natl Acad Sci U S A 114: 11721-11726. 869 Turner JM. 2007. Meiotic sex chromosome inactivation. 870 Van der Auwera GA, O'Connor BD. 2020. Genomics in the cloud: using Docker, GATK, and 871 WDL in Terra. O'Reilly Media, Inc. 872 Vibranovski MD, Zhang Y, Long M. 2009. General gene movement off the X chromosome in the 873 Drosophila genus. Genome Res 19: 897-903. 874 Vicoso B, Charlesworth B. 2006. Evolution on the X chromosome: unusual patterns and 875 processes. Nat Rev Genet 7: 645-653. 876 Vicoso B, Emerson JJ, Zektser Y, Mahajan S, Bachtrog D. 2013. Comparative sex chromosome 877 genomics in snakes: differentiation, evolutionary strata, and lack of global dosage 878 compensation. PLoS Biol 11: e1001643. 879 Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, Cuomo CA, Zeng Q, 880 Wortman J, Young SK et al. 2014. Pilon: An Integrated Tool for Comprehensive 881 Microbial Variant Detection and Genome Assembly Improvement. PLOS ONE 9: 882 e112963. 883 Wang W, Jia M, Zhao C, Yu Z, Song H, Qin Y, Zhao W. 2021. RNF39 mediates K48-linked 884 ubiquitination of DDX3X and inhibits RLR-dependent antiviral immunity. Science 885 advances 7: eabe5877. 886 Wang Z, Zhang J, Yang W, An N, Zhang P, Zhang G, Zhou Q. 2014. Temporal genomic 887 evolution of bird sex chromosomes. BMC Evol Biol 14: 250. 888 Webster TH, Vannan A, Pinto BJ, Denbrock G, Morales M, Dolby GA, Fiddes IT, DeNardo DF, 889 Wilson MA. 2024. Lack of Dosage Balance and Incomplete Dosage Compensation in the 890 ZZ/ZW Gila Monster (Heloderma suspectum) Revealed by De Novo Genome Assembly. 891 Genome Biol Evol 16: evae018. 892 W estfall AK, Telemeco RS, Grizante MB, Waits DS, Clark AD, Simpson DY, Klabacka RL, 893 Sullivan AP, Perry GH, Sears MW et al. 2021. A chromosome-level genome assembly 894 for the eastern fence lizard (Sceloporus undulatus), a reptile model for physiological and 895 evolutionary ecology. Gigascience 10: giab066. 896 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint 31 Wu CI, Xu EY. 2003. Sexual antagonism and X inactivation--the SAXI hypothesis. Trends Genet 897 19: 243-247. 898 Xu L, Irestedt M, Zhou Q. 2020. Sequence Transpositions Restore Genes on the Highly 899 Degenerated W Chromosomes of Songbirds. Genes (Basel) 11: 1267. 900 Yang Y, Chang T-C, Yasue H, Bharti AK, Retzel EF, Liu W-S. 2011. ZNF280BY and ZNF280AY: 901 autosome derived Y-chromosome gene families in Bovidae. BMC Genomics 12: 13. 902 Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Molecular biology and 903 evolution 24: 1586-1591. 904 Yonenaga-Yassuda Y, Rodrigues MT, Pellegrino KCM. 2005. Chromosomal banding patterns in 905 the eyelid-less microteiid lizard radiation: The X1X1X2X2:X1X2Y sex chromosome 906 system in Calyptommatus and the karyotypes of Psilophthalmus and Tretioscincus 907 (Squamata, Gymnophthalmidae). Genet Mol Biol 28: 700-709. 908 Yu Y-H, Lin Y-W, Yu J-F, Schempp W, Yen PH. 2008. Evolution of the DAZ gene and the AZFc 909 region on primate Y chromosomes. BMC Evol Biol 8: 96. 910 Zhou Q, Zhang J, Bachtrog D, An N, Huang Q, Jarvis ED, Gilbert MTP, Zhang G. 2014. 911 Complex evolutionary trajectories of sex chromosomes across bird taxa. Science 346: 912 1246338. 913 Zhu Z-X, Matsubara K, Shams F, Dobry J, Wapstra E, Gamble T, Sarre SD, Georges A, Graves 914 JAM, Zhou Q et al. 2022. Diversity of reptile sex chromosome evolution revealed by 915 cytogenetic and linked-read sequencing. Zool Res 43: 719-733. 916 917 .CC-BY-NC-ND 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 7, 2025. ; https://doi.org/10.1101/2025.04.01.646716doi: bioRxiv preprint