Analysis of coverage in M. picta and M. parae males and females, and search for genes with diploid coverage in M. picta males

Microsatellite genotyping indicated hemizygosity of chromosome 12 for all loci tested, in both M. picta and M. parae [18]. Hemizygosity was confirmed for M. picta by coverage data obtained by high-throughput sequencing approach (SeqSNP, see Methods). The full details are in a S1 Methods file. The SeqSNP data yielded an estimate of 99% gene loss from the Y [18].

Our new analyses of coverage, used whole genome low coverage sequences of multiple males and females sampled from multiple natural populations. The Methods sub-sections “Estimating divergence between Poecilia and Micropoecilia sequences” and “Analysis of coverage in M. picta males and females”, describe our reference-guided local assembly approach to map M. picta sequences to genes in the P. reticulata female reference genome sequence [5]. The analysis included many more genes than SeqSNP, and confirmed that most guppy LG12 genes show haploid coverage in M. picta males, while autosomal genes show diploid coverage levels. Fig 2 shows the normalised coverage male/female (M/F) ratios in the two sexes for chromosome 12 and, for comparison, an autosome (LG15). The ratios for the autosomal genes were indeed consistently close to 1 for the vast majority of the 824 LG15 genes where this could be estimated (the median value was 1.01 and only 5 genes yielded values below 0.7). Other autosomes gave very similar results (see Dryad accession number https://doi.org/10.5061/dryad.cjsxksn65); the results are similar to those already published for M. picta [19]. For LG12, the ratios are close to 0.5, indicating hemizygosity, for 822 genes whose orthologue positions in the guppy assembly span most of the chromosome (median 0.53). There appear to be slightly more outliers for chromosome 12 than for the autosome, possibly reflecting transposable element or other repetitive sequences; however, as the sequences analysed were genes, there should be few such cases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Coverage depths of genes on the M. picta sex chromosome and an autosome.

In each part, the y axis shows estimated normalized coverage values and the x axis shows physical positions in the guppy reference assembly of the homologous chromosome. In both the guppy and M. picta, chromosome 12 is the sex chromosome, and LG15 is shown as a representative autosome (as explained in the text, other autosomes are similar). The red horizontal line indicates a value of 0.5, as expected for hemizygous genes, and the blue one indicates diploid coverage (a value of 1), as seen for most genes across the autosome, and for two clusters of LG12 genes.

https://doi.org/10.1371/journal.pgen.1009704.g002

Male hemizygosity for LG12 genes appears to be species-wide in both M. picta, including in our M. picta males from Trinidad samples where M. parae is absent, as well as in all M. parae populations tested [19,24,25]. We found the same state in a third Micropoecilia species, M. bifurca, based on genotypes of six informative microsatellite markers spanning most of the sex chromosome (S6 Table).

We next searched the M. picta LG12 gene coverage results for any genes that might be retained on the Y chromosome of this species, using the same natural population samples. Genetic degeneration is expected to produce consistent signals of hemizygosity across all populations of a species, and lack of degeneration should consistently yield similar depths of coverage in both sexes (see Methods and S3A Fig). Sequences with male-female coverage ratios close to 1 could represent either PAR genes in M. picta, or genuinely non-hemizygous parts of the fully sex-linked region (so called “XY” sequences), or autosomal sequences that have been mis-assembled on the guppy LG12. This will affect our assemblies, as the guppy genome assemblies were used to assign positions for our assembly of the M. picta sequences.

In contrast to the majority of LG12 genes described above, 47 of 57 genes with positions distal to 25 Mb in the guppy female assembly (from 25,067,709 to 26,240,873 bp) had M/F coverage ratios close to 1 (the median value for these 57 genes = 1.0); their terminal locations in the guppy genome assemblies suggest that these might be PAR genes in both species. In addition, 16 of the 21 genes with positions from 9,954,185 to 10,927,843 bp in the female guppy assembly, also had high values (the median for this set was 0.95). The values for these sets both differed significantly from the ratios of genes in the hemizygous regions, though 12 genes at the terminal end of the guppy assembly (26,262,815 to 26, 423,538 bp) had lower values (median 0.53). These may represent an error in the guppy assembly, as their genetic map positions are not consistent with the positions assigned in the female guppy assembly [11], and the positions have been changed in the male assembly, based on PacBio sequencing long reads, which also does not support an inverted order of the X or the Y [6]. It seems likely that these may not to be PAR genes, though definitive evidence will require genetic mapping with denser genetic markers. Given their locations and diploid coverage values in males, the set of 47 terminal region genes seem likely to be PAR genes in M. picta. In the next section, we describe genetic mapping results supporting this conclusion.

Before doing so, we evaluate the set of 16 genes that are assigned positions near 9–10 Mb in the guppy female assembly. As noted above, these might be autosomal genes that have been mis-assembled, or further PAR genes, rather than fully Y-linked in M. picta, but not hemizygous (so called “XY” sequences). Several observations support the mis-assembly hypothesis. First, these genes’ positions in the guppy female assembly are not well supported. This region is repetitive, and is much more AT-rich than its neighbouring regions in the guppy female assembly (S4 Fig), and the gene sequences have low mapping quality. Second, consistent with assembly problems, these sequences are placed in a different location in the guppy male assembly, where they appear in a centromere-proximal region, between 2. 6 and 3.1 Mb [6], but they map near the end of the platyfish (Xiphophorus maculatus) assembly that corresponds with the terminal, pseudo-autosomal, part of the guppy assembly (S5 Fig). Comparison with the platyfish assembly finds a set of sequences from this LG12 region (including several genes, and a repetitive sequence) on Xm1, the homologue of the guppy autosome, LG7 (S6 Fig). An LG7 location is also supported by the observation that a single long contig in our M. picta whole genome sequence [26] carries these genes along with homologues of many guppy LG7 genes. Overall, therefore, our detailed coverage analyses did not detect any good candidates for M. picta genes with retained Y-linked alleles.

Unlike the hemizygous genes, which gave consistent signs of hemizygosity in most populations, genes that appeared to be duplicated were mostly detected only in samples from individual populations (S3B Fig). At most, 32 genes suggested consistent evidence from all populations, and these were not studied further. We did, however, perform a further test for genes in the terminal part of LG12 that might have Y-linked or other duplications, using BLAT searches of M. picta contigs with all genes in the guppy PAR and adjacent regions as the queries, and excluding targets that were not single-copy in M. picta. The few gaps in the region distal to 20 Mb tend to correspond with regions of low gene density, and to gaps between contigs, and do not indicate duplicated genes in the most terminal region, with the highest M/F ratios (S7 Fig).

Moreover, we note that the coverage ratios in the putative PAR are not all due to coverage in males exceeding that for autosomal genes, as Fig 2 might suggest, but frequently reflect lower coverage in females, relative to autosomal genes. S8 Fig shows these results for the part of the M. picta LG12 corresponding to the region shown in S7 Fig, relative to the values for an autosome, LG15. Although we excluded low coverage and low-quality sequences (as described in S1 Methods 1), genes in the region between 20 and 24 Mb (where the M/F coverage ratio remains close to 0.5), have low relative (M/A and F/A) values than for the rest of LG12. This suggests possible problems with mapping reads in this region. In the 64 genes in the most distal region, where M/F coverage values are often, only ten genes have higher relative coverage in males than females, 14 have similar coverage in both sexes, and 40 have lower values in males, and the overall median M/F ratio is 0.98. These results suggest highly variable success in mapping reads from PAR sequences, perhaps due to their high GC content [17].

Genetic tests of sex linkage in the guppy and M. picta

The observation that coverage is similar in males and females for most M. picta homologues of guppy PAR genes suggests that the M. picta sex chromosome pair has not evolved into a X0 system, but retains a small PAR region carrying genes homologous to those in the guppy PAR. Alternatively, this region might no longer be part of the M. picta sex chromosome, i.e. these genes might be autosomal in this species.

To test whether these genes are sex-linked or autosomal in this species, and in M. parae, we studied two genes from the terminal region of the guppy sex chromosome, zer1 and fpgs (see Methods); in the guppy genome assembly, both are distal to the most distal fully sex-linked marker genetically mapped [3]. We next describe results demonstrating that both genes show sex-linkage, and must on be located on LG12 in all three species. To distinguish between complete and partial sex-linkage, we also genotyped samples from natural populations.

Linkage analysis of the zer1 gene in the guppy QHPG5family (see Methods), based on length differences of intron 10 of this gene, confirms that this gene is partially sex-linked in the guppy. It recombines with the sex-determining locus in male meiosis, and its genetic map location is between two other pseudoautosomal markers CTT98 and GT354 (S7 Table of ref. [18]), also included in the Dryad files for the present manuscript). The fpgs gene did not have markers suitable for mapping in the guppy.

In M. parae, variants in families suggest sex linkage of the fpgs gene, though the family sizes are small. The zer1 gene was not polymorphic in this species and therefore could not be mapped (S5 Table of ref. [18]). Genotypes in natural populations, however, suggest that both these genes appear to be partially sex-linked, as in the guppy (S4 and S5 Table). Among 24 individuals genotyped for fpgs intron 4 length variants, seven of the 13 males were heterozygous, and two alleles were shared with the female sample. Three males were heterozygous for an allele not seen among the females, but only 5 females were present in the sample, so these are probably merely low frequency alleles, not male-specific. Similarly, many males were heterozygous for length variants of the zer1 gene intron 10 shared with females (in total, 5 different alleles were detected in this small sample of males, two of which, plus one other allele, were found in the small female sample).

In M. picta, the zer1 gene is sex-linked in families from Trinidad and Suriname rivers, based on variable intron 10 sizes. In all families whose sires were not homozygous for the allele with length 957-bp, only the male progeny inherited the father’s other allele (families T1, T2, T6, T7 and S3 in S4 Table). However, genotypes in samples from natural populations suggest complete sex linkage (S5 Table). All 31 Trinidad females were homozygous for a 957 bp allele of zer1 (also present in the Suriname sample, but not fixed there). However, among 31 Trinidad males, only 18 were homozygous for this allele; the other 13 were heterozygous for one of at least eight male-specific alleles, suggesting that zer1 does not recombine with this species’ sex-determining locus in male meiosis. The 957-bp allele thus appears to be carried by all X chromosomes, while the Y is highly polymorphic for many different Y-specific alleles. SNPs in the fpgs gene also show sex-linkage in male meiosis (family S3 in S4 Table). In the Trinidad natural population samples, all males appeared to be homozygotes for the microsatellite marker in this gene, including two males with rare variants, one of which was also present in a female heterozygote, suggesting that males are hemizygous, and that this gene is also completely sex linked in M. picta (S5 Table). The broods from all seven Trinidad females also all appeared to be sired by homozygotes; however, two Suriname sires were inferred from their progeny genotypes to be heterozygous for the SNP marker in the gene, and therefore this region cannot be hemizygous (S4 Table). Overall, the data suggest that the zer1 gene, and possibly the more distal fpgs gene, are completely sex-linked in M. picta, but partially sex-linked in the guppy and M. parae (Table 1). Some apparently hemizygous males might simply reflect heterozygtes for common null alleles for the fpgs gene microsatellite. The M. picta Y chromosome may therefore have stopped crossing over between the XY pair in a region that still recombines in the other two species.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Results of genetic mapping or genotyping of two guppy PAR genes in two Micropoecilia species.

https://doi.org/10.1371/journal.pgen.1009704.t001

Sequence divergence between Poecilia and Micropoecilia species and the time when genetic degeneration occurred in Micropoecilia

Finally, we wished to test between the two possible evolutionary origins of the guppy Y chromosome outlined in the Introduction section above. Because no genes seem likely to be Y-linked copies that have escaped loss during the degeneration in M. picta Y chromosome, we cannot compare the time of Y-X divergence with the time since the split of M. picta and M. parae from the lineage leading to P. reticulata. However, we have estimated synonymous site divergence (K_s values, see Fig 3) between M. picta and P. reticulata for 27 genes, from several different chromosomes (see Methods). The median K_s value is 0.055. If the M. picta Y stopped recombining and started to degenerate after the split of the two lineages (possibility 1 of Fig 1), its complete degeneration, and evolution of equal expression of sex-linked genes in males and females, indicating dosage compensation [19], must have occurred in an evolutionary time corresponding to at most half this value, or a K_s value of about 2.8%. Given that our microsatellite genotyping suggests that the M. parae and M. bifurca Ys are also highly degenerated (see results above and ref. [24]), and that synonymous site divergence between the Micropoecilia species indicates a split at least half as long ago as the split between the Poecilia and Micropoecilia lineages (as illustrated in Fig 3), it is likely that the degeneration occurred in the Micropoecilia lineage before it split from the Poecilia lineage, rather than after the Micropoecilia species split. This conclusion is conservative, as the only other possibility if the ancestor had a non-degenerated Y (possibility 1 above) is that degeneration occurred in one Micropoecilia species, and the degenerated Y introgressed into the other species. This implies complete degeneration during a time when the three species were able to hybridise, most likely early in the time after the split from the Poecilia lineage. In the next section, we assess whether the time is long enough for complete genetic degeneration to evolve.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. DNA sequence divergence between the guppy and M. picta and between two Micropoecilia species.

Divergence estimates are for synonymous sites in coding regions of 27 genes, with Jukes-Cantor correction for saturation of substitution, and 95% confidence intervals of the estimates. The horizontal lines indicate the median values (5.4% for P. reticulata versus M. picta, and 4.1% for M. picta versus M. parae. Divergence estimates between all three Micropoecilia species are given in S3B Table.

https://doi.org/10.1371/journal.pgen.1009704.g003

How degenerated is the M. picta Y-linked region?

The evidence described above confirms that the M. picta Y-linked region is highly degenerated. Our analysis suggest that the Y carries at most 60 genes, compared with about 900 on the X, and that most of these are orthologous to genes at the terminal (PAR) end of the guppy chromosome 12 in either the female or male assembly.

It is not yet completely excluded that the M. picta Y chromosome might carry sequences similar to their X-linked sequences, with the Y-linked ones too diverged to be detected by previous analyses [7,18,19]. Highly diverged Y-linked gene sequences could fail to map to the relevant location on the P. reticulata reference X. To maximise our ability to detect any Y-linked sequences that are still present, our analyses, high-throughput genotyping with primers designed from P. reticulata [18], and mapping to genome sequences, described here, targeted coding sequences, which are least likely to diverge. The mean divergence of M. picta sequences from guppy sequences for all site types, which is the most relevant for the success or failure of PCR, is only 2.7% (see S3 Table), and indeed most targeted genes in M. picta were successfully amplified in our high-throughput genotyping experiment. Any undetected gene sequences retained on the Y of this species must be more diverged than this average.

Our coverage analysis using genome sequences of males should have detected some of these diverged sequences, if they indeed exist. The read mapping employed was very sensitive, but was quite intolerant of mismatches, and would therefore be unlikely to call Y-linked sequences that are extremely diverged. Genes whose X cognates are missing from the reference X (notably in difficult to assemble, repetitive regions) would also remain undetected. However, despite employing a P. reticulata reference sequence, our mapping had a high overall recall value at high quality (for the 13 females from the CAR_H population alone, we recovered 140 out of an initial original 180.5 million reads, or 78%, see S1 Methods 1). This conservative statistic suggests that only small amount of raw sequence was too diverged to map. The larger number of genomes used for the coverage analysis included much higher numbers of reads, and recall was similar for both sexes. A better assembly of M. picta in the future should indicate whether some regions have so far escaped inclusion in assemblies, and might yield an even higher recall value, and further reduce the proportion of unmapped sequence. Our detailed analyses of LG12 detect fewer than 100 genes, many of which are probably PAR genes. This is consistent with the inferred number of guppy PAR genes, based on genetic mapping of genes with known physical positions in the assemblies [3], or on high GC content of their introns [17]. As described above, PAR genes are also present in M. parae. Neither Micropoecilia species appears to have an X0 system.

Under the proposed scenario in which the guppy Y evolved from an X chromosome (possibility 2), the boundary of the highly recombining PAR1 of chromosome 12 might differ from the PAR boundary in Micropoecilia species. The evidence from the zer1 gene (see Results) suggests that this boundary might indeed differ, and that the M. picta PAR could be slightly smaller than the guppy one. Further data are needed to determine the PAR boundaries in both species, and to test whether any difference reflects a change in M. picta.

Overall, we conclude that (apart from PAR genes) the M. picta Y has probably retained at most a small proportion of coding sequences ancestrally present on the X. To estimate the extent of gene loss from the M. picta Y since it stopped recombining with the X chromosome, we can assume that the guppy and M. picta X chromosomes carry similar numbers of genes, as our results show that the M. picta PAR does not include many more genes than the terminal guppy PAR1. On this basis, the X-linked region carried about 800 ancestral genes when recombination stopped in the Micropoecilia lineage, of which at most 16, or 2%, have been retained on the Y.

Both our results [18] and genome sequencing [24] suggest that Y chromosome degeneration is also advanced in M. parae. The finding of male hemizygosity in all three Micropoecilia species tested (see Results and S6 Table), supports the suggestion that degeneration probably evolved in a common ancestor before these species split [24]. The variant version of possibility 1 (outlined above), involving introgression of a degenerated Y into two other species, is therefore not likely.

More specifically, recombination probably ceased in a common ancestor, though degeneration might have proceeded independently in the three Micropoecilia lineages. As microsatellite sequences are in non-coding regions and introns, the hemizygosity of LG12 microsatellites in males suggests that parts of this chromosome have been deleted and that only a remnant of the ancestral Y remains. This conclusion remains to be tested by cytogenetic observations of the physical size of the M. picta and M. parae Ys.

How long has the M. picta Y been non-recombining, and is it long enough for complete genetic degeneration?

The lack of any good candidate genes with both X- and Y-linked copies described above prevents our estimating Y-X sequence divergence to understand whether the M. picta Y stopped recombining with the X chromosome before the split from the Poecilia lineage, and, if so, how long before. However, it is still possible to evaluate the two possibilities for the origin of the guppy sex chromosomes summarised in Fig 1, because, under possibility 1, degeneration in the Micropoecila species started after their split from Poecilia lineages (Fig 1A). We can therefore assess whether there was enough time for almost complete degeneration to occur. The K_s value of 2.7% mentioned above over-estimates the time of the split, because it is a “raw divergence” value that includes an unknown contribution from polymorphisms within the ancestral species [27]. The synonymous site diversity in the ancestor is unknown, and there are as yet no estimates from the species of interest, so we cannot currently correct for such polymorphism and estimate the net divergence value that estimates the separation time in terms of fixed substitutions since the time of the split [27]. However, even conservatively using the raw divergence value, it seems unlikely to represent enough time for the documented almost complete gene loss from the fully Y-linked region, and evolution of dosage compensation [19]. Given that all three Micropoecilia species so far studied appear to have degenerated Y chromosomes, it is unlikely that the Y chromosome stopped recombining in a common ancestor and degeneration occurred independently in each lineage after their split. Based on synonymous site divergence estimates (S3B Table), this restricts the available time to about half that suggested by the time of the guppy-M. picta split.

The evolutionary time needed for such changes is not yet fully understood, but we next outline current predictions and empirical information from animals (data from plants are less helpful, because selection in the haploid pollen may cause degeneration to be slow, as reviewed in ref. [28]).

Empirical data are very scarce on the relationship between the extent of degeneration, and divergence at synonymous sites (which, assuming a rough molecular clock, relates well to the number of generations since recombination stopped). The best data are from humans. However, all regions of the human Y have become highly degenerated, with the Y-linked genes either being completely lost, or present but non-functional [29], including the two young strata in humans that are pseudo-autosomal in other mammals, including a lemur species [30,31]. These two strata are not very young, as they are present in lemur species [32], and have mean Y-X K_s values of 45% and 14% (the estimates are rough, because these strata include only 17 X-linked genes [29]). The human strata are already too old to provide information about partial degeneration over times corresponding to that for the Micropoecilia lineage degeneration.

The currently available data of which we are aware for degeneration over smaller evolutionary times are listed in S7 Table. Although data are scanty, all species with Y-X (or W-Z) K_s > 20% show complete degeneration, but degeneration is incomplete when K_s < 20%, including the very well-studied D. miranda neo-Y chromosome [33]. Data are available for a few fish fully sex-linked regions, and these are incompletely degenerated, despite having K_s considerably higher than the value corresponding to the time available in Micropoecilia. The oldest stratum in the threespine stickleback, Gasterosteus aculeatus (with about 600 genes, slightly fewer than in the species analysed here), has an estimated X-Y K_s = 16% (or a time back to the Y-X split corresponding to K_s = 8%, assuming the same neutral substitution rate for both X- and Y-linked regions) yet 18% of X-linked genes still have detectable Y-linked copies. This fully sex-linked region is shared with G. wheatlandi, but in that species only 8% of genes are found on the Y, and these yielded a lower K_s estimate, 9% [34]. In both species, many Y-linked genes are still present in regions that stopped recombining more recently; in G. aculeatus, regions with Y-X K_s = 3–4% (corresponding to an evolutionary time similar to that available for degeneration of the Micropoecilia Y) have retained around 70% of genes [35]. Even during double this time, gene loss might conservatively be expected to be at most 60% (double the 30% value estimated for the stickleback), given that initial rapid degeneration quickly slows down as the number of remaining functional genes undergoing deleterious mutations declines [23]. In an unrelated fish, C. semilaevis, 44% of fully sex-linked genes are present on the W as well as the Z, indicating much less degeneration than in M. picta, but the W-Z divergence time is considerably longer, with W-Z K_s estimated to be 15% [36]. The time available therefore seems unlikely to be compatible with the almost complete gene loss from the M. picta Y, especially if the degeneration occurred before the split from other Micropoecilia species, as argued above.

Theoretical modelling provides another source of information. However, the time needed for complete degeneration of a Y chromosome cannot yet be predicted, and many uncertainties remain. One set of parameter values modelled led to 90% degeneration of a Y-linked region with 2,000 genes in about 5 million generations, but with < 1,000 genes degeneration is greatly slowed (compare the right-hand panels of Figs 8 and 4 of ref. [23]).The K_s value of 0.054 between the guppy and M. picta estimated above suggests about 27 million generations since the Micopoecilia species split from the Poecilia lineage (assuming a neutral mutation rate of 10⁻⁹ per site per generation), which might be enough time for nearly complete degeneration, given about 800 ancestral genes (see above). Degeneration can be faster with smaller population sizes, higher deleterious mutation rates, stronger selection against deleterious mutations, or unequal fitness effects of strongly deleterious mutations [23], but allowing mutations on the X in a diploid population can slow it down [37]. The genetic degeneration process initially involves Muller’s ratchet combined with background selection, both of which occur fastest when the effective population size is small. A population bottleneck could therefore explain faster degeneration in a species such as M. picta, compared with other species. Neutral diversity has not yet been estimated in any of the fish species studied here, and effective population size estimates are not yet available. Importantly, however, this part of the process slows down after about half of the genes are lost from a fully Y-linked region, and cannot account for complete degeneration (see Fig 3 of Bachtrog 2008). The same figure illustrates that genetic hitch-hiking as advantageous mutations spread in a fully sex-linked region is slower in small populations, though this (or its combination with the ratchet, as shown in Fig 8 of the same paper) can eventually lead to complete degeneration. Assuming no back mutations (which would slow down the ratchet process), 90% loss of genes occurs in about 5 million generations, in a population with 10,000 Y chromosomes, and about 98% after 15 or 20 million generations. Complete degeneration would take longer in a larger population size.

Overall, the limited data suggest that there have probably been too few generations since the split between the Poecilia and Micropoecilia lineages for the extensive degeneration observed. However, more quantitative empirical data on degeneration and further modelling work are needed to better understand the time required for degeneration to become complete.

Could the guppy Y have evolved recently from an X chromosome?

The evolution of the guppy Y from an X chromosome (as proposed in the second possibility above) is not implausible. First, turnovers in sex-determining genes are well documented in several groups of fish. Among species related to those discussed here, the platyfish sex-determining locus is on the homologue of the guppy chromosome 20 [38], and the more distantly related genus Oryzias is well known for turnovers [39], which have also been documented in cichlid fish [40], which are more distantly related still [41]. Second, the housefly offers a precedent for movement of a male-determining factor onto the X chromosome [42,43]. Mobile sex-determining regions are also known from other organisms, including Salmonid fish [44].

Turnovers are, however, thought to be less likely to occur when a species has a highly degenerated Y or W chromosome, because they may lead to low fitness homozygotes for these chromosomes being generated [45]. The proposed change in the lineage ancestral to the guppy might thus seem implausible. However, a turnover that involves a new sex-determining factor on the same chromosome, as proposed here, does not generate such homozygotes. Any slight advantage to the new system can therefore lead to a turnover [45]. A possible advantage (in an XY system, such as that in the ancestor proposed here) is that loss of genes from a degenerated Y may reduce the fitness of heterozygous males, so that the evolution of a new sex-determining factor on the X might be favoured. This is similar to the advantage proposed for repeated turnovers driven by degeneration after sex-linked regions appear [46–48].

Loss of the ancestral Y chromosome would clearly be prevented if it carried genes essential for male fertility or viability, and this has been suggested as an explanation for the long-term evolutionary stability of sex chromosomes in groups such as mammals and birds [45]. As the M. picta fully Y-linked region appears to carry very few genes, it may not carry such male function genes. Another characteristic that might make a system like that in M. picta unlikely to undergo a turnover event is that it has evolved dosage compensation of sex-linked genes [19]. If a male-determining factor appears on an autosome or an X chromosome and creates a new male heterogametic system, with loss of the ancestral Y, males in the new population will have two X chromosomes. If the dosage compensation mechanism involves higher expression of X-linked genes in males, these “XX males” might experience a deleterious high expression levels of these genes (potentially higher than in females, which will presumably have evolved to close to optimal levels). However, the mechanism controlling dosage compensation need not respond to the phenotypic sex. If instead, it responds to X chromosome dosage, the expression level could be the same in the new male genotype as in females (and in the ancestral XY genotype with hyper-expression from the X). Although M. picta has clearly not evolved an X0 sex-determining system, with the X-autosome balance controlling development of the individual as phenotypically male or female, dosage compensation might be controlled in this way.

A final potential barrier to the turnover proposed here is that the conditions for a turnover become more restrictive if sexually antagonistic genes have accumulated on the ancestral Y chromosomes [49]. However, as these authors point out, evolution of sex-limited expression of such genes on the ancestral sex chromosomes reduces any sexual conflict and allows SA alleles to become fixed, rather than remaining polymorphic, increasing the probability that a new sex-determining factor in a different genome region can invade. The absence of Y-linked coloration polymorphisms in M. picta, unlike the guppy and M. parae, suggests that this barrier may also not apply to a system like that in M. picta, though clearly the ancestral state is unknown with respect to this characteristic. Overall, the evidence from M. picta suggests that several factors expected to inhibit turnovers may not apply.

How did the guppy and M. picta Ys evolve?

If the ancestor of the guppy and M. picta had a Y that was already highly degenerated before the split between Poecilia and Micropoecilia, a duplication of a pre-existing autosomal or sex-linked gene could have been involved, or a new maleness factor could have arisen [50]. Both types of turnover are known in fish (reviewed in refs. [39,51]). Unlike these cases, which were detected because the XY pairs are not homologous in the different species, the guppy sex chromosome pair has not changed, implying that an X-linked gene mutated, or arose by a duplication from another genome region, as in the medaka, O. latipes. In the guppy, such a duplication could have originated on an autosome or the ancestral X or Y. For example, an ancestral Y-linked male-determining factor might be superseded by a mutation of a gene that was present on the ancestral X, that conferred a male-determining action. If the ancestral Y did not carry a male-determining factor, and the ancestor had a balance system, a mutation in an X-linked gene is more plausible than duplication of or mutation in a Y-linked gene. A duplication might have led to suppressed expression of an X-linked gene that represses male functions in “normal” females, and thus created a male phenotype. Duplications with phenotypic effects that arise by repressing the expression of the progenitor gene are documented in several organisms, and is termed “co-suppression”, first discovered in transgenic plants and animals [52,53]. A similar phenomenon occurs in sex-determining systems of plants in the family Salicaceae [54–56].

It is not known what mechanism stopped crossing over between the Y and X of M. picta and M. parae, but we speculate that, if chiasmata in the ancestral species were localised to the chromosome termini, as in the guppy [3], even a small Y chromosome inversion might have prevented crossing over between the XY pair, as diagrammed in Fig 4. An inversion could explain the complete complete sex linkage in M. picta of some genes in the guppy LG12 terminal region, PAR1 (as the rearrangement’s terminal boundary would determine the PAR boundary in stage 2 of Fig 4). However, loss of recombination in this region probably occurred recently, after recombination was suppressed across most of LG12 in M. picta, as two genes are in a non-recombining region in this species, but not in the guppy or M. parae (Table 1). Furthermore, Y-X divergence (based on intron sequences in the zer1 gene whose sizes identify Y-versus X-linked sequences), is only 0.6%, so the PAR boundary probably changed recently in M. picta, similarly to changes detected in mammals. Finally, if a sexually dimorphic recombination pattern (heterochiasmy) was already established in a species ancestral to M. picta and P. reticulata, it would offer a simple explanation of how the ancestral X chromosome could have evolved into a Y in stage 4 in Fig 4 (possibility 1 described in the Introduction section).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Diagram of hypothesis 2 for the evolution of the guppy Y chromosome.

This hypothesis involves the appearance of a new maleness (M) factor on an X chromosome like that in M. picta. The diagram shows 4 proposed stages, indicated by numbers, from a hypothetical ancestral state (numbered 1) to the guppy state (4). Regions where crossovers occur are indicated by ‘x’ symbols, with high and low densities of symbols indicating high or low crossover rates, respectively. The ancestral state might rarely have undergone crossovers, indicated by ‘x’ symbols in a small distal part of the centromere-proximal end of the chromosome pair. A completely non-recombining Y-linked region, like that in extant M. picta males, could have evolved by an inversion (stage 2, see text), followed by complete genetic degeneration (stage 3). Possible sources of a new maleness factor in the guppy are discussed in the text. If phenotypic males rarely have crossovers in the part of the chromosome to the left of the maleness factor, this region would become closely linked with maleness. In the case of a duplication, a completely Y-linked region might also be created (expanding the guppy Y, compared with the X), but the centromere-proximal part of the chromosome could still occasionally recombine with the X, preventing genetic degeneration.

https://doi.org/10.1371/journal.pgen.1009704.g004

The proposed appearance of a maleness factor on the X during the evolution of the guppy Y resembles events proposed in Muscid flies [42]. In some Dipteran flies, the maleness factor can appear on different chromosomes [57], including the X [43,58]. In Dipterans, crossing over does not occur in male meiosis, so that such a maleness factor immediately causes the recipient chromosome to be transmitted without recombination with the X. In the guppy, crossover patterns are sexually dimorphic, with crossovers strongly localised to the chromosome termini in males [3]. If male meiosis in the ancestor was similar, as our results suggest is true in the close outgroup of the species studied here, the platyfish, Xiphophorus maculatus [17], a new Y could arise in a manner similar to that in the house fly. However, Muscid fly sex chromosome turnovers appear to involve a fusion with a physically small ancestral chromosome (Muller element F) that carried an ancient sex-determining factor and is missing as an independent chromosome from such species, rather than appearance of a new maleness factor [42]. The guppy sex-determining factor is currently unknown, but a fusion scenario seems unlikely, because the male genome assembly [6] provides no evidence that any part of the XY pair is derived by fusion with, or transposition from, another chromosome, and no large region is present on the guppy XY pair but absent from either the M. picta X [26] or from the guppy Y [6]. It was recently suggested that the guppy sex-determining region arose by a duplication of a chromosome 9 region [59], but a careful examination did not support this, and suggested that assembly error could be involved [11].

It is also unlikely that the guppy Y-linked male-determining factor could have arisen by an X-Y crossover event. Such events occasionally occur in male guppies [11,15,60]. However, if the ancestral Y was degenerated, like the M. picta or M. parae Y, this implies complete lack of recombination, as a small recombination rate prevents degeneration [46].

A final possibility is that a mutation or duplication occurred on the ancestral X. Small duplications have been detected in the guppy sex chromosome pair, including in the LG12 region assembled near 20 Mb [6,11]. However, these this might simply be sequences that are linked to the male-determining locus, and not the proposed new maleness factor.

Ethics statement

The University of Edinburgh School of Biological Sciences ethics committee, SBS Ethics, reviewed and approved this research (Ethics Assessment Reference: dcharles-0001).

Fish samples, DNA extraction

M. picta samples were collected by L Yong, at Cunipia River, part of the Caroni swamp, Trinidad (10° 36 ’ 18.89 “N, 61° 25 ’ 28.462 “W), and in several other Trinidad rivers (S1 Table) by Cameron Ghalambor. A mixed natural population of M. picta and M. parae was sampled at Pronkweg road, Commewijne District, Suriname (GPS coordinates N 0712221, E 0659682). Males of the two species can be distinguished by their coloration phenotypes with reasonable accuracy [71], but females are difficult to distinguish. We therefore used several molecular markers which were found to be diagnostic for the two species (the microsatellite AG179, and intron size variants of the pseudo-autosomal genes fpgs and zer1 (see S4 and S5 Tables). Genomic DNA for microsatellite genotyping, genotyping of intron length variants, and for high-throughput genotyping (see below), was extracted using the Echolution Tissue DNA Kit (BioEcho, Germany). M. bifurca samples from Suriname were provided by Dr. David Reznick. The identities of all Micropoecilia samples were verified by sequencing their cytochrome oxidase genes, which show synonymous site divergence of at least 12% (and the sequences obtained all exceeded 437 bp).

The samples were collected under permit number 150/19 from the Ministry of Agriculture, Animal Husbandry and Fisheries of Suriname, and under an un-numbered permit from the Fisheries Division of the Department of Agriculture, Land and Fisheries of the government of Trinidad and Tobago. The study protocols were reviewed and approved by the University of Edinburgh School of Biological Sciences (see https://apps.bio.ed.ac.uk/ethrev/#/view/dcharles/1).

Microsatellite genotyping (using primers listed in S2 Table), together with initial inference of gene loss in M. picta and M. parae are described in Bergero et al. [18], which describes coverage data obtained by a high-throughput genotyping (SeqSNP) experiment. The experiment included 10 M. picta individuals of each sex, and was carried out by LGC Genomics (LGC Genomics GmbH, Ostendstrale 25, 12459 Berlin, Germany, www.lgcgroup.com/genomics); it yielded counts of the numbers of sequences obtained for each targeted region. The regions targeted were selected from within coding sequences, with about 50 bp of sequence flanking each region also being coding sequence, in order to maximise the chance that the sequence would amplify in diverse populations, and to minimise the representation of repetitive sequences. The SNPs and their locations in the guppy genome assembly are listed in a file deposited in Dryad (see above). Sequences were included only if all individuals yielded a genotype, and the mean count for the targeted site exceeded 10 per individual (usually greatly in excess of this threshold); one targeted sequence was excluded because of an exceedingly high count, suggesting a repetitive sequence. Overall, 1,338/1,749 (77%) of the autosomal coding regions targetted yielded sequences satisfying our criteria that all M. picta individuals yielded an allele count, and 411/491 of the LG12 regions (84%). The high success rate supports the expectation that nucleotide divergence between M. picta and the guppy is low (see next section and the Results section). New analyses of coverage, based on complete genome sequences, are described below.

Estimating divergence between Poecilia and Micropoecilia sequences

As no reliably assembled genome is available, to estimate sequence divergence between M. picta and other species, we used whole genome sequences obtained from Cameron Ghalambor as raw read data (see S1 Methods 1). The names of all populations, and the sample sizes are in S1 Table. The methods used for analysing these genome sequences are described in detail in the S1 Methods text.

To obtain coding sequences from these data, an extended bioinformatics protocol was enacted using a reference-guided local assembly approach similar in principle to software such as SRAssembler and the SHORE/SUPERLOCAS low-coverage resequencing pipeline (S1 Methods 1). Because of the low depth of coverage, we used de novo assembled reads that mapped to genes in the P. reticulata reference genome sequence, as previous analyses of a small number of nuclear genes had shown that sequence divergence between M. picta and P. reticulata is not high [72]; for synonymous sites, the sequences available in GenBank yielded a divergence value of about 7%, and divergence based on all site types of about 2.3% (S3B Table). S1 Methods 2 describes the methods used for obtaining the coding sequences.

The between-species divergence estimates used data from M. picta females from a population sampled from a high salinity (10–19 PPT) region of the Caroni river, Trinidad (population name CAR_H). Coverage for each female was low, with a mean of 13.9 million reads, or 1.4x based on the size of the P. reticulata genome, but the pooled data from the thirteen females with high coverage and reliable sexing represented 180.6 million reads, or 18.5x coverage. The alignments for each individual gene were analysed using MEGA v7 [73], using the complete deletion option. We estimated synonymous and non-synonymous divergence, and divergence for all site types, between the species included (K_s, K_a and K_{all_sites} and the sequence lengths are given in S3 Table). The MEGA files will are available in Dryad.

Analysis of coverage in M. picta males and females

The low coverage genome sequences of our hundreds of M. picta samples were from individuals whose genomes were sequenced were of unknown sexes. We therefore used the genome sequence data to determine the sexes of the sampled individuals, based on the hemizygosity of most M. picta X-linked genes in males (S1 Methods 3, see also S2 Fig). S1 Methods 4 describes how the coverage data were used to search for M. picta LG12 genes with consistently diploid male/female coverage ratios, and to test for genes that might have additional Y-linked copies, as indicated by M/F ratios > 1.4. We also examined the coverage results for possible duplications onto the M. picta Y (S3A Fig and S3B Fig).

Genetic and population genetic analyses in P. reticulata and in Micropoecilia species

We used two approaches to test for partial sex-linkage of markers that are assembled in the guppy sex chromosome terminal region (which are not hemizygous in in M. picta males, see the main text), genetic mapping and genotyping individuals sampled from natural populations.

First, we mapped two genes in a P. reticulata family whose parents were from a high-predation locality in the Quare river., family QHPG5 of ref. [3]. Both parents of this family were genotyped previously, and a genetic map was estimated for several microsatellite markers and numerous SNPs, distributed across LG12. For genetic mapping, we developed markers in two terminal region genes, zer1 (at position 26,208,835 bp, encoding a zer1 homologue protein), and fpgs (at position 26,286,204, which encodes a folylpolyglutamate synthase protein, mitochondrial isoform X2). The total estimated physical length of the guppy chromosome 12 is 26.5 Mb, and both these genes are distal to the most distal fully sex-linked marker in the guppy genetic maps estimated from the QHPG5 family and other families [3,11]. The zer1 gene was mapped using a length variant in intron 10, and fpgs using SNPs and length difference variants for intron 4.

The zer1 and fpgs genes, along with microsatellite markers used to determine individuals’ sexes [18], were also genotyped in samples of fish from natural populations, and in families, from M. picta and M. parae. The numbers of genes analysed was limited because several genes and microsatellite markers from the 1-Mb terminal region of the guppy sex chromosome assembly failed to amplify in M. picta, or yielded multiple bands, and could not be analysed.