Polygenic sex determination in the cichlid fish Astatotilapia burtoni

Comparative scaffolding of the A. burtoni genome

A. burtoni was one of five cichlid species sequenced as part of the cichlid genome project, but was also the species with the most fragmented genome assembly (scaffold N50 of 1.2 Mb, versus 2.5–4.4 Mb for the other four species; [41]). In order to provide chromosome-level organization of A. burtoni scaffolds for mapping, we ordered and oriented them to anchored genomes of two of other cichlid species (Additional files 1 and 2). The Nile tilapia (Oreochromis niloticus) genome assembly was anchored to a high-resolution radiation hybrid map, providing a golden path [41, 42]. Genome scaffolds of Metriaclima zebra, a Lake Malawi species, were also anchored to a RADseq-based genetic map [43]. A. burtoni scaffold assignment was consistent by linkage group for both species, but regular inconsistencies in scaffold order indicate a moderate level of intrachromosomal rearrangement between O. niloticus and M. zebra. We use the O. niloticus scaffolding for the following analysis because it provided the most parsimonious ordering based on our marker data, including producing continuous, uninterrupted mapping intervals. We number A. burtoni linkage groups based on synteny with established O. niloticus linkage group numbering, per convention in the cichlid genetics community [42].

Identification of two sex determination loci

Two A. burtoni families with near-equal sex ratios were genotyped using a ddRADseq strategy, and resulting SNP marker genotypes were tested for association with gonadal sex (Table 1). Markers with significant p values (p < 0.01) were further sorted by allele frequency in each sex, to search for minor alleles at intermediate frequency in one sex, and absent in the other, following a dominant model of sex determination (Additional file 3). Both families had numerous markers with significant associations with phenotypic sex, and alignment of marker sequences to the tilapia genome revealed clustering of the associated markers by linkage group. In one family, marker genotypes in a region of LG13 show near perfect association with sex, as do markers on unanchored portions of the tilapia assembly including UNK25 and UNK107 (Fig. 1, Additional file 3). Based on association and haplotype patterns we predict that sex-associated scaffolds aligning to UNK25 and UNK107 should be placed on LG13, within or adjacent to the sex associated scaffolds aligning there. Because maternal alleles at LG13 co-segregate with sex, we designate this family as having sex determined by a ZW system.

Family	M	F	Readsa	Stacksa	Polymorphic stacksb	Lowest p-valuec	Corresponding Tilapia linkage groupsd
RAD 1	12	12	9.7 × 106	435,006	4,777	9.6 × 10−6	LG13, UNK25, UNK107
RAD 2	9	8	5.9 × 106	341,384	4,268	1.4 × 10−3	LG5, LG14

Families with number of male (M) and female (F) offspring. ^aAverage reads and stacks per offspring. ^bTotal number of polymorphic stacks used for association mapping, requiring a minimum stack depth of 10 to call homozygotes, and genotypes called for 20 of 24 offspring for RAD 1, and 15 of 17 offspring for RAD 2. ^cThe lowest p-value obtained for Fisher’s exact test of allele association with sex, shared by multiple markers (Additional file 3). ^dLinkage groups with markers with lowest p-value association

In a second family, associated RAD markers consistently aligned to tilapia LG5 and LG14 (Fig. 2a, Additional file 3). Previous karyotype analysis demonstrated that A. burtoni has two fused chromosome pairs relative to all other haplochromine cichlids surveyed [44]. Examination of our associated marker genotypes at LG5 and LG14 revealed clear haplotypes with no recombinants spanning the two linkage groups (Additional file 4), demonstrating that LG5 and LG14 constitute one of the chromosome fusions in the A. burtoni karyotype relative to tilapia; we refer to this fusion as LG5-14 in A. burtoni (Fig. 2). Paternal alleles co-segregate with sex at LG5-14, and thus we designate it as carrying an XY sex determination system.

In families where LG13 is associated with sex determination, either a maternal or a paternal allele can be found associated with sex, depending on the family tested (Table 2). In several families, a maternal haplotype segregates with sex in the offspring, and the same maternal haplotype is associated with female sex in multiple families, supporting the presence of a dominant W female sex determination allele (Additional file 6). Surprisingly, paternal alleles at LG13 co-segregate with sex of offspring in families from two sires, suggesting that a dominant Y male sex determination allele may also be present on LG13 (Table 2, Additional file 6). Differences in parental inheritance of sex determination alleles among families with LG13 sex determination suggest that an XYW sex determination system akin to that found in platyfish may be present at LG13; however, occasional sex reversal could also account for observed patterns in the absence of a LG13 Y allele (see below).

Several families showed no association of sex with markers at LG13 or LG5-14 (Table 2). These families allow the possibility of additional, unidentified sex determination loci in A. burtoni. Alternatively, interactions among identified sex determination alleles, or between sex determination alleles and the environment, could confound interpretation.

Interactions between sex loci

Because we find multiple loci and alleles that appear to impact sex, we would expect some families to be segregating multiple alleles that would interact to influence sex determination. Studies in both tilapias and Lake Malawi cichlids with polygenic sex determination have revealed clear patterns of epistasis relating to phenotypic sex [21–23, 30]. Unfortunately, we find no clear examples of such interactions in our families. One small family (Family 19) has perfect segregation of sex with paternal alleles at LG13, supporting the sire having an LG13 XY genotype. Two sons but no daughters in the family also inherit a multi-marker haplotype at LG5-14 associated with the LG5-14 Y in other families (Additional file 7). We hypothesize that the genotype of the sire and two sons in this family is LG5-14 XY, LG13 XY, and thus that individuals with a Y allele at each locus are viable and develop as males; however, this interpretation is based on a single, small family and is thus very tentative.

In families tested where sex is not clearly associated with alleles at either LG13 or LG5-14 alone (Table 2), inter- and intra-locus interactions may be present that do not conform to our expected patterns of allelic dominance. Additionally, some families were not included in this study because they had extreme sex ratios (sometimes single-sex), and/or few offspring, preventing meaningful interpretation of inheritance patterns; however, interaction of multiple sex alleles may produce skewed sex ratios or non-viable genotypes. In cichlids and other fish species, co-segregation of both Y and W alleles in families produces female-biased sex ratios [18, 22, 30], and in pygmy mice with an XYW system, mating of some sex genotypes results in inviability of a portion of their offspring [45]. Parnell and Streelman described five distinct sex determination loci (three XY and two ZW systems) in a Lake Malawi cichlid mapping cross, demonstrating that sex ratios represent the outcome of complex interactions between relatively strong and weak sex determination alleles [23]. Thus, epistasis may be difficult to characterize, particularly in small families, and further work is needed to understand the phenotypic outcomes of interactions between the sex determination alleles identified.

Sex ratios and penetrance

Previous studies have demonstrated that polygenic sex determination systems in cichlids produce skewed sex ratios. In a survey of species from Lake Malawi, families with sex determined by an XY system had equal sex ratios, while those with a ZW system generally had a female-biased sex ratio [22]. In families where both XY and ZW systems segregated, a 1:3 (male:female) sex ratio is expected [22]. In some Lake Victoria cichlids, sex determination factors on supernumerary chromosomes also produce female-skewed sex ratios [32]. Comparing the sex ratios of our LG5-14 XY families and LG13 ZW A. burtoni families reveals an intriguing and significant pattern: XY families have a male-biased sex ratio, and ZW families have a female-biased sex ratio (Fig. 3a). Because we only identified a few putative LG13 XY families, we do not test for significant divergence in sex ratio or penetrance.

Comparing penetrance of each sex determination allele, we find no significant differences (Fig. 3b), and a mean penetrance of 88 % for all sex determination alleles. Incomplete penetrance could be explained by a number of factors, including incorrect phenotype data, recombination between our sex-linked markers and the sex determination polymorphism, or sex reversal. If sex reversal is the case, it presents an alternative hypothesis regarding the interpretation that both W and Y alleles are present at LG13. The presence of a W allele has greater support based on replication with a number of dams, but only two sires were found to have their LG13 alleles segregate with sex of their offspring. Phenotypic males in some families with LG13 ZW dams occasionally carry a W haplotype at tested markers. For example, see Family 3, where two males exhibit the W-associated haplotype at LG13 (Additional file 6). Uncommon males that are genetically female (ZW at LG13) would have their LG13 alleles co-segregate with the sex of their offspring, lending the appearance of an XY system even though a distinct, male-determining Y allele might not exist at the locus.

Genetic intervals for sex loci

Incomplete penetrance makes it difficult to determine mapping intervals via recombination breakpoint mapping, because it is unclear if genotype-phenotype mismatch at a particular marker is due to recombination or sex reversal. Though our ddRADseq families included relatively few individuals, we estimate genetic intervals based on strongly associated genome scaffolds. For the LG13 ZW locus, scaffolds spanning approximately 4.6 Mb show the highest association with sex (Fig. 1, Additional file 3). At LG5-14, a broad plateau of association spans 43 Mb of the linkage group based on data from the initial discovery family (Fig. 2). We further interrogated the association data with a sliding window analysis on LG5-14, to look for regions with the highest density of associated SNPs, which may correspond to the region of highest sequence divergence. Sliding window analysis reveals a broad peak of the highest density of associated SNPs on several scaffolds that align to LG5 of tilapia. In one family with LG5-14 XY sex determination, recombination breakpoints in one individual support that the region under this peak contains the sex determination locus (Individual 8, Family 18, Additional file 6; spanning scaffold 156 to scaffold 66).

Candidate genes

Though precise genetic intervals are not clear, we examined gene content on the A. burtoni scaffolds with the highest association in our RAD family data (Additional file 8). Notably, genes with established roles in sex determination in multiple vertebrate species (e.g., sox3, dmrt, amh, and gsdf) are not found on these scaffolds, suggesting the use of novel sex determination genes by A. burtoni. Here we highlight intriguing candidate genes on associated scaffolds based on known roles in sexual development or fertility.

A second chromosome fusion

Because our ddRADseq haplotype data was able to resolve the linkage groups involved in one of two the A. burtoni chromosome fusions relative to other cichlids [44], we attempted to identify the other known fusion by exploring linkage between markers in our RAD genotype data. Building an outcross linkage map with only 24 offspring from RAD family 1 produced a poorly resolved linkage map, with 30 linkage groups containing greater than five markers. Though the map was very rough, it identified candidates for linkage group fusions that were then examined by manual comparison of haplotypes. Markers aligning to two tilapia linkage groups, LG8-24 and LG16-21, were linked in our map output. (The nomenclature for both LG8-24 and LG16-21 in tilapia does not indicate fused chromosomes, but is instead a historical artifact: the standard O. niloticus linkage map originally had 24 linkage groups [73], but later improvement of the genetic map identified 22 linkage groups [42], collapsing the former LG8 and LG24 into a single linkage group, as well as LG16 and LG21.) We find continuous haplotypes spanning LG8-24 and LG16-21, with significant linkage between markers aligning to the two linkage groups (LOD = 10.3, Additional file 4), supporting the fusion of the two associated chromosomes in A. burtoni relative to other cichlids. Maintaining the community naming convention based on the original O. niloticus linkage map, this linkage group would be given the rather cumbersome designation LG8-24-16-21. A previous report suggested that one of the fusions in A. burtoni consisted of LG15 and LG19 relative to other cichlids [31], but we find no evidence of linkage between markers mapping to those linkage groups (LOD = 0.8, Additional file 4). This discrepancy may stem from interchromosomal rearrangements more complex than simple fusions, where different portions of a particular tilapia LG map to different A. burtoni LGs. Thus, use of different genetic markers could provide patterns suggesting fusion of different LGs.

Comparisons among cichlid sex determination systems

A previous study inferred the mode of sex determination system in A. burtoni using hormonal sex reversal [40], in a line derived from the same lab strain we use here. Estradiol sex reversal of males and subsequent testcrosses strongly indicated an XY sex determination system, and ruled out a ZW system [40]. These results are consistent with our mapping data, because it is quite possible that only the XY sex determination system was present in the four families used for sex reversal experiments. Indeed, bottlenecks during transfer of the A. burtoni lab strain between labs may inadvertently lead to complete loss of some sex determination alleles. Concurrent mapping experiments by another group using the same lab line as the sex reversal experiments above identified co-segregation of sex with markers on LG5 indicating an XY system [74], completely consistent with the mapping results we describe here. While the concurrent study found no association of sex with LG13, it did identify an association of sex with LG18 in the offspring of a pair of wild-caught individuals from a Lake Tanganyika estuary [74]. Together, our differing but overlapping results support the presence of complex, polygenic, and potentially rapidly evolving sex determination in A. burtoni.

A number of sex determination loci have been identified in cichlids. In tilapia, sex determination has been associated with LG3, LG20, LG23, and two distinct loci on LG1, depending on the species and population [30, 35, 36, 75]. In Lake Malawi haplochromine cichlid species, LG3, LG5, LG20, and two loci on LG7 have been associated with sex determination [21, 23]. In haplochromine cichlids from Lake Victoria, LG2, LG5, and supernumerary B chromosomes have been associated with sex [31, 32]. Here we identify two loci impacting sex determination in the riverine haplochromine A. burtoni, LG5-14 and LG13. While this is the first time LG13 and LG14 have been associated with sex determination, it is intriguing that LG5 has now been associated with haplochromine sex determination in A. burtoni, Malawi, and Victoria cichlids. To provide finer scale comparison of sex association with LG5, we identified A. burtoni scaffolds containing sex-associated markers from other studies (Fig. 2b), revealing overlap of sex loci among multiple species. As on-going mapping studies identify specific polymorphism underlying sex determination in each context, we should be able to distinguish reuse of ancestral polymorphism versus repeated de novo evolution of sex determination on LG5. However, both W [21, 22] and Y (this study; also suggestive patterns in two families in [22]) acting alleles have been associated with LG5, suggesting at least two evolutionary events. LG5 has been dubbed a “speciation chromosome” because of the number of adaptive and species-delineating phenotypes mapped there, including pigmentation, vision, jaw shape, and tooth shape [76], and it may exhibit dynamic evolution of sex determination with both W and Y alleles now clearly identified. The identification of the LG5-14 fusion in A. burtoni should provide additional context to consider the role of this chromosome in adaptation and speciation, particularly in comparing the evolutionary trajectory of LG14 among lineages with and without the fusion. Indeed, the LG5-14 fusion in A. burtoni may represent karyotype evolution driven in part by the presence of a sex determination locus on LG5. For example, if novel alleles under strong sexually antagonistic selection arose on LG14, the resulting sexual conflict could have been resolved by fusion to the LG5 sex chromosomes, ensuring that sex-specific fitness benefits were only expressed in the appropriate sex.

Evolutionary context

Evolutionarily, the genus Astatotilapia represents an immediate outgroup to the adaptive radiations of cichlid species in Lake Malawi and Lake Victoria [41, 77], providing important context for studies of sex determination evolution in those flocks. Ongoing studies may reveal if the sex determination locus on LG5/LG5-14 is the first found both within and without the lake adaptive radiations, and/or shared among radiations. Additionally, a recent report implicated the same region of LG5 as one of a handful of genomic islands of speciation in Astatotilapia species, based on sequence divergence between ecomorphs ([78]; Fig. 2b). One hypothesis is that the signal of divergence at LG5 in that study is due to divergence related to the sex determination locus, with the LG5 XY system potentially present in one ecomorph and absent in the other.

Ultimately, our mapping results must be confirmed in natural populations. The recent discovery that the natural sex determination chromosome in zebrafish (Danio) was lost in lab strains urges caution in interpreting genetic mapping of sex determination in domesticated fish [28]. More recently, large-scale changes in gene expression were identified after only a single generation of domestication in trout [79], suggesting that domestication can act rapidly on the phenotypic output of the fish genome. The standard lab strain of A. burtoni we use here has been maintained in the lab for nearly four decades [80], representing over 100 generations in captivity in which loss or gain of sex determination alleles could have taken place. Indeed, recently collected stocks of A. burtoni exhibit a number of phenotypic differences from the classic lab strain [81], suggesting underlying genetic differences between wild and lab strain A. burtoni.

In general, the evolutionary pressures that maintain polygenic sex determination systems remain unknown. We favor the hypothesis that the multiple sex genotypes produced by polygenic sex determination may provide alternative fitness strategies, for example, alternate sex genotypes may be associated with differences in color (yellow vs. blue male morphs) or behavior (potential sneaker vs. territorial mating strategies) in A. burtoni [82, 83]. In any case, A. burtoni should serve as an excellent model for understanding the outcomes of polygenic sex determination for primary and secondary sexual characteristics. The availability of genome editing tools in the species permits mechanistic dissection of genetic variation [84, 85], while decades of research in the species provides a wealth of information regarding sexually dimorphic behavior, physiology, and related gene expression [38, 81, 86–89].