Structure and decay of a proto-Y region in Tilapia, Oreochromis niloticus

Sequencing of male and female DNA pools

We obtained ~202 million reads from the pool of male DNA and ~219 million reads from the pool of female DNA. 90.12% of the male and 90.67% of the female reads were aligned to the O. niloticus reference genome. Genome-wide coverage was slightly lower in males (32.97, standard deviation = 24.41 alignments per site), compared to females (36.68, standard deviation = 31.39 alignments per site).

Large block of divergence on linkage group 1

The mean F_ST between the male- and female-pooled genomes at polymorphic sites over the entire genome was 0.0356 (standard deviation = 0.030). A region between 10.1 Mb and 18.9 Mb on linkage group 1 showed a substantially higher value of F_ST = 0.0807 (standard deviation = 0.061) (Figures 1a, and 2a). This region corresponds to the previously mapped sex-determination region in this strain of O. niloticus[21, 24–27]. Mean read coverage within the differentiated region was lower in males (34.65, standard deviation = 10.56), compared to females (38.45, standard deviation = 12.00), but this difference was consistent with the total number of reads obtained from each sex. We used Fisher’s exact test to determine whether the allele frequency of SNPs was significantly different between males and females. We found a cluster of highly significant SNPs within the differentiated block on linkage group 1 (Figures 1b, and 2b).

We also counted the number of positions per 10 kb window that were fixed in female pools and had a SNP in intermediate frequency in male pools, as would be consistent with females having two X chromosomes and males having an X and a Y chromosome, using Sex_SNP_finder_now.pl. There were 40,514 of these SNPs found across the genome. 18,277 (2,076.932/Mb) lay inside the differentiated block and 22,237 (24.197/Mb) lay outside. Among the 300 non-overlapping 10 kb windows with the highest frequency of these SNPs, 290 were found within the differentiated block on linkage group 1. The mean number of such SNPs per window was 21.81 (standard deviation = 13.84) within the differentiated block and only 0.33 (standard deviation = 1.29) outside of this region (Figures 1c, and 2c). The elevated F_ST, along with the abundance of intermediate frequency SNPs in males that are fixed in females, suggests that this region has limited, if any, recombination between the X and Y alleles.

We considered the possibility that this block of differentiation is an artifact of the process by which we selected individuals for sequencing. We initially screened individuals by genotyping two sex-linked microsatellites in order to confirm family identity and sex. We required males to demonstrate a heterogametic pattern and females to demonstrate a homogametic pattern for both markers. Five male and five female individuals were excluded by these criteria and may represent naturally sex-reversed individuals. The sharply defined edges of the block lie 4.22 Mb upstream and 3.37 Mb downstream of the microsatellites we genotyped (Figures 1 and 2), which would normally represent approximately 5 cM of genetic distance in this species [28]. However, there is no evidence of an exponential decay of F_ST in the flanking regions as would be expected if there was recombination between the markers and the edges of the block. We also considered the possibility that the high level of differentiation might be due to an 8.8 Mb duplication on the Y. However, the depth of read coverage is relatively consistent across this entire linkage group. Additionally, cytogenetic studies have not revealed any evidence of heteromorphy in this chromosome pair as would arise from a translocation [29]. The sum of the evidence suggests that this block of differentiation most likely reflects an 8.8 Mb inversion on the Y-chromosome.

The relatively small size of the putative inversion, and its location in the middle of the chromosome, make it challenging to characterize using standard cytogenetic techniques. Ideally, we would characterize the breakpoints, but we were unable to identify anomalous Illumina mate pairs near the ends of the inversion in our short insert libraries. Longer reads or more widely spaced mate pairs will be needed to characterize the breakpoints of the proposed inversion.

Functionally significant SNPs

The elevated density of high impact SNPs within the proposed inversion leads us to believe that deleterious alleles have begun to accumulate on this proto-Y. This is in accordance with the canonical model of heterogametic sex-chromosome evolution [2, 32] and empirical observations of the therian mammal Y-chromosome, Silene, Drosophila and tongue sole [4, 33–35].

Localization of the sex-determining gene

Previous studies have concluded that sex is multifactorial in O. niloticus[24, 36] with a major sex-determination gene on LG1 [21]. Our study confirms this previous work, identifying an XY sex-determination locus in the middle of LG1 (Figure 3). The sex-determination gene was first mapped near microsatellite markers GM201 (13.96 Mb) and UNH995 (18.02 Mb, Figure 3a) [24]. Additional AFLP and FISH mapping found sex-associated markers at 13.79 Mb, near 18 Mb and at 19.43 Mb (Figure 3b) [25, 26]. Another study confirmed GM201 and UNH995 along with several other sex-associated markers spanning a region from 7.05 Mb to 18.02 Mb (Figure 3c) [21]. Lastly, a RAD-seq experiment found the highest associations at 14.95 Mb (LOD score 18.5), but demonstrated a broad region spanning 10.92 Mb to 16.44 Mb with a LOD score above 15 (Figure 3d) [27].

The multifactorial nature of sex-determination in this species causes difficulties for genetic mapping studies. An XX individual may develop as a male due to other genetic factors, or environmental effects on differentiation. These individuals would appear to be recombinant in the sex interval. We previously claimed to exclude Wilm’s tumor protein homolog (Wt1b) as the sex-determining gene on the basis of two recombinant individuals [37], but this conclusion is now in doubt. Conversely, the absence of recombination within the proposed inversion may preclude any further fine-mapping of the gene responsible for sex-determination.

Differences in gene expression

The block of differentiation on linkage group 1 comprises just more than 1% of the assembled genome and contains 257 RefSeq annotated genes. Cufflinks predicted 234 gene models within the block of differentiation and predicted 22,411 gene models across the entire transcriptome. Of the gene models that showed an FPKM of >0.05 in at least one sex, 7,977 gene models (37.4%) showed higher expression in males, while 13,375 (59.7%) gene models showed higher expression in females. Furthermore, within the inverted region, only 68 of these gene models (29.6%) showed a male bias (Additional file 3), while 162 of these gene models (69.2%) showed a female bias (Additional file 4). The enrichment of female biased gene models within the proposed inversion is significant (χ² = 5.58, p <0.05). These data suggest that this sex chromosome is at an early-to-intermediate evolutionary phase where the degradation of a proto-Y has begun and expression of Y-linked genes in males is reduced. However, mechanisms for complete dosage compensation have yet to take hold.

Candidate sex determiners

First, we analyzed all SNPs that SnpEff classified as high impact mutations (Table 1). One prominent candidate within the proposed inversion is Ras-related protein R-Ras2 (10.49 Mb-10.51 Mb), which is part of the Ras-MEK-ERK pathway within the TGF-β signaling network [38]. Alterations to the TGF-β network have been suggested as the mechanism for sex-determination in several fish species [15]. Ras2 has been implicated as particularly important in the proliferation of cells [39] and is expressed during early development in a hermaphroditic fish, Kryptolebias marmoratus[40]. Ras-related protein R-ras2 has a stop codon gain in intermediate frequency in males that is absent in females. Disruption of R-ras2 could lead to decreased cell proliferation of primordial germ cells, resulting in increased likelihood of maleness [15, 39].

Next, we evaluated SNPs that SnpEff categorized as missense mutations (Additional file 2). The first of these candidate genes is Wilms tumor protein homolog, Wt1b (14.86 Mb-14.88 Mb), which has been implicated in gonadal development and acts directly upstream of AMH, the sex-determination gene in Odontesthes hatcheri[41]. Wt1b has also been demonstrated to bind to DNA and upregulate the sex-determination gene Sry in mammals. There is an A237V missense mutation in Wt1b that is absent in females and in intermediate frequency in males. Although our previous paper rejected Wt1b on the basis of two recombinant individuals [37], in light of the proposed inversion, we now believe that these individuals represented instances of natural sex reversal, not recombination.

A third candidate gene is estrogen-related receptor gamma, ERRγ (17.05-17.11 Mb). It has a R172H missense mutation within a predicted DNA-binding domain [42]. ERRγ has been shown to be a transcriptional activator of DAX-1, and DAX-1 has been implicated as having an antagonistic effect to Sry in mammals [43]. Therefore, a mutation in the DNA-binding domain of ERRγ could reduce DAX-1 transcription and thus have a masculinizing effect.

Growth regulation by estrogen in breast cancer 1 (GREB1) is another candidate gene (17.41-17.42 Mb) with a missense mutation. The R1775C mutation alters the side chain from a basic side chain to a polar side chain. GREB1 has been shown to be predominantly expressed within ovaries of young mice [44]. Additionally, GREB1 has been demonstrated to be a coactivator of estrogen receptor-α [45]. Therefore, the missense mutation in GREB1 could downregulate the expression of estrogen receptor-α, resulting in a masculinizing effect on the developing embryo.

Another potential sex-determination gene is transcription factor SOX-6 (10.22 Mb-10.30 Mb). There is a T789K missense mutation in intermediate frequency in males that is fixed in females and changes a polar side chain into a basic one. SOX-6 protein is localized to the same nuclear speckles as Sry and it has been suggested that it might play a role in sex-specific splicing in mammals [46].

We also evaluated gene models showing differential expression between males and females (Additional files 3 and 4). AFG3(ATPase Family Gene 3)-like protein 1 (13.72 Mb-13.73 Mb) has over a nine-fold male-biased expression. It is also on the list of SNPs with high impact coding alterations with a stop codon gain. However, a clear tie to sex-determination has yet to be elucidated.

Suppression of tumorigenicity 5 protein (10.80 Mb-10.85 Mb) and Ras association domain-containing protein 10 (11.40-11.41 Mb) were also identified for having over a three-fold male-biased expression pattern. Ras association domain family proteins have been implicated as tumor suppressors [47–49]. Therefore, upregulation of these genes could suppress primordial germ cell proliferation leading to maleness.

Lastly, it is possible that there could be Y-specific genes that were not captured in our study, because the reference genome that the reads were aligned to is a homozygous clonal XX individual.

Genomic DNA pools

All animal procedures were conducted in accordance with University of Maryland IACUC Protocol #R-10-73. The fish sequenced are 3rd generation descendants of fish collected from a commercial tilapia farm in Amherst, Massachusetts USA. Individuals from two related lab-raised families were sacrificed and visually inspected for the presence testes or ovaries to determine the sex of each fish. Fish with ambiguous or immature gonads were excluded from the study. DNA was extracted from fin clips using a standard phenol/chloroform protocol. To confirm the family identity we genotyped each individual for two sex-linked microsatellite markers selected from the Broad anchored tilapia assembly on linkage group 1 (MS1045 at 14.32 Mb and MS1141 at 15.53 Mb). 33 males and 20 females from family BYL078 and 25 males and 13 females from family BYL084 were selected for pooling. DNA from each individual was then quantified by Picogreen fluorescence on a BioTek FLx800 spectrophotometer and appropriate dilutions were made to ensure equal representation of each individual in the pooled samples. The pooled male (or female) DNA from each family was sheared to 500 bp using a Covaris shearer and indexed separately during library construction. Paired-end libraries for each family/sex were constructed for Illumina sequencing using the Illumina TruSeq DNA Sample Prep Kit (Illumina, San Diego, CA). The male (or female) libraries from each family were combined and each sex was sequenced in a separate lane on an Illumina HiSeq 2000. The male and female reads were deposited to NCBI with the accession numbers SRR1606298 and SRR1606304, respectively. Only reads passing the Illumina CASAVA filtering were retained.

Read qualities was checked with FASTQC [50]. Alignments to the O. niloticus anchored reference assembly [51] were performed with Bowtie 2 [52] using the --very-sensitive setting (Additional file 5). The mean alignment rate was 90.12% in males and 90.67% in females (Separate values for each family are given in Additional file 5). Read alignments were filtered for a minimum mapping quality (MAPQ) of 20 before further analysis. Insert sizes were analyzed using Picard tools CollectInsertSizeMetrics package [53]. The aligned mean insert size was 188.76 bp (standard deviation = 44.81 bp) for males and 167.62 bp (standard deviation = 37.29 bp) for females. Variants were called using GATK [54].

Genomic analysis

Popoolation2 [55] was used to calculate F_ST and Fisher’s exact test on allele frequency differences between the male and female pools. Initial F_ST results from the individually adapter-indexed families were very similar, so all subsequent analyses were performed on the combined male or female pool, including unassigned reads from the male and female lane which could not be assigned to a particular family.

A custom Perl script, Sex_SNP_finder_now.pl (available at https://github.com/Gammerdinger/sex-SNP-finder), was used to identify SNPs at intermediate frequencies in the male pools and were fixed or nearly fixed in female pools at the same position. Intermediate SNPs were defined as SNPs with a frequency between 0.3 and 0.7 within the male pool. Fixed or nearly fixed sites required a frequency less than or equal to 0.1 or greater than or equal to 0.9 within the female pool. We used a non-overlapping window of 10 kb to determine the density of these SNPs. The non-overlapping window did not include positions with coverage less than 10 reads in both sexes. The Sex_SNP_finder_now.pl script outputs a tab-delimited file with the number of SNPs per window along with an Integrative Genomics Viewer file [56] that lists all SNPs that were fixed or nearly fixed in one designated pool and in intermediate frequency in the other.

We used SnpEff [30] to identify variants predicted to alter gene function. The SnpEff output was filtered to consider only the SNPs found using Sex_SNP_finder_now.pl. SnpSift [31] was used to extract out SNPs with similar effects and impacts. A complete list of genes within the proposed inversion can be found in Additional file 6.

Transcriptome analysis

Gonads were dissected from individual larvae 28 days post-fertilization. The sex of each larvae was determined by genotyping microsatellite markers highly associated with sex. RNA from approximately 20 male or 20 female larvae was pooled and cDNA libraries were constructed using the Illumina TruSeq DNA Sample Prep Kit. Sequencing of these libraries yielded ~392 million reads for each male and female pool. Reads were aligned to the O. niloticus reference sequence with TopHat2 [57]. NCBI RefSeq annotations were used to guide the Cufflinks [58] assembly (-g) and Cuffdiff was used to was used to determine FPKM values for those gene models. The results were subsequently filtered to exclude gene models whose FPKM value was less than 0.05 in both males and females. Additionally, when comparisons between FPKM of the two sexes was carried out, if the FPKM value exceeded 0.05 in one sex and was zero in the other sex, it was considered an undefined bias favoring the sex with expression. Female-biased and male-biased gene models from inside and outside the proposed inversion were counted and statistical significance was looked for using χ² with Yates’ correction on a 2x2 contingency table. These male and female reads from the RNA-Seq experiment were deposited to NCBI with the accession numbers SRR1606274 and SRR1606273, respectively.