BACKGROUND
the differentiation of an evolving pair of sex chromosomes involves many evolutionary steps, including changes in gene content on the x and y  <cit> ; gene loss and the fixation of deleterious mutations on the y  <cit> ; epigenetic modifications related to silencing of y-linked genes  <cit>  and the evolution of dosage compensation  <cit> .
the drosophila miranda neo-sex chromosome system is emerging as a useful resource to study the evolution of young sex chromosomes  <cit> .
d. miranda split from its sister species d. pseudoobscura only about  <dig> my  ago, and has sex chromosomes of different ages: muller a , the ancestral sex chromosome shared among all drosophila species including d. melanogaster, is at least  <dig> my old  <cit> ; muller d , a neo-sex chromosome system of intermediate age  and shared with d. pseudoobscura; muller c , specific to the d. miranda lineage.
the genome sequence of d. miranda has recently been published, including a draft assembly of the neo-y chromosome  <cit> .
a major challenge in the assembly of this chromosome was the fact that the neo-y has already accumulated many repetitive sequences, which are difficult to assemble – similar to the ancestral y chromosome of drosophila which has not been assembled at all.
in addition, the high level of sequence identity between the neo-x and neo-y especially in coding regions further complicates assemblies using next-generation sequencing data, since many short sequencing reads are identical from the neo-x and neo-y.
to annotate the neo-y and study expression changes, diagnostic snps between the neo-x and neo-y were used to estimate the relative abundance of neo-sex linked transcripts  <cit> .
this, however, does not allow studying the structure of the neo-y transcripts per se.
in particular, it is unknown how complete neo-y protein-coding genes are, especially with regards to the transcripts that contain premature stop codons or frameshift mutations.
furthermore, the use of diagnostic snps for expression analysis will lead to an inherent mapping bias towards higher neo-x expression, since rna-seq reads can overlap with parts of transcripts that differ structurally between the neo-x and neo-y.
using the neo-x sequence as a baseline for comparison, one can only study neo-y transcripts that are also present on the neo-x and show no major structural changes or indels.
thus, there is considerable interest to study the neo-sex transcriptome directly, which we do here.
in addition to providing insights into sex chromosome evolution, the neo-sex chromosomes of d. miranda also allow us to study other cellular phenomena, such as nonsense-mediated mrna decay .
nmd is a cellular surveillance mechanism that degrades transcripts containing premature termination codons   <cit> .
destruction of transcripts containing ptcs by nmd prevents production of truncated, potentially harmful proteins that may interfere with normal cellular processes.
the exact mechanism through which nmd recognized ptcs in drosophila is unclear, though the length of the 3′ utr has been implicated  <cit> , whereas the presence of an exon-exon junction downstream of a stop codon is an initiator of nmd in mammals  <cit> .
given that a large number of neo-y transcripts contain ptcs, the d. miranda system makes an excellent case for testing if the length of the 3′utr leads to nmd in flies.
a major challenge in constructing the neo-sex transcriptome is that coding regions on the neo-x and neo-y are very similar; based on the genome sequence, the estimated divergence at coding regions is only about  <dig> %.
accordingly, using illumina rna-seq reads for transcriptome assembly, not all reads can be assigned to the neo-x or neo-y unambiguously.
this problem is for the most part irrelevant for the neo-x transcriptome, which can be assembled directly in females .
in males, however, transcriptome assemblers fail to assemble the neo-x and neo-y into separate transcripts.
instead, chimeric transcripts are produced that contain partial neo-x and partial neo-y sequence, connected by regions that are not differentiated between the neo-sex chromosomes.
to obtain an assembly of the neo-y transcriptome, we developed a bioinformatic pipeline, making use of genomic read mapping against the neo-sex transcripts in males.
methods
assembly pipeline
we used illumina genome analyzer ii paired-end 100-bp rna-seq reads of the inbred line msh <dig> to assemble the d. miranda transcriptome with trinity  <cit> .
first, low-quality sections of rna-seq reads and remaining adaptor sequence were discarded, and the minimum contig size of the trinity assembly was set to 200 bp.
for the neo-y assembly,  <dig> million rna-seq reads from four different samples  were combined in a single assembly; similarly, to obtain neo-x transcripts,  <dig> million rna-seq reads from virgin ovary, virgin female whole body and virgin female gonadectomized carcass were combined and assembled together .
for the autosomes and old sex chromosomes , both male and female reads were combined for a single assembly.
the d. pseudoobscura transcriptome  was assembled de novo using combined virgin male and virgin female whole body samples .
an overview of the pipeline to assemble the neo-y transcripts is given in figure  <dig>  using blastn, trinity transcripts in males were compared to the d. miranda genome assembly  <cit> .
potential muller c transcripts  were kept and modified to become neo-y transcripts, by introducing variants that distinguish the neo-x from the neo-y.
for this means, male genomic reads were compared to the potential muller c transcripts and snp positions were identified.
note that, in contrast to rna-seq reads, male genomic reads contain equal proportions of neo-x and neo-y variants in males.
to avoid introducing erroneous variants at intron-exon junctions, male trinity transcripts were first divided up into exons.
this was done by mapping male genomic reads against the male muller c transcripts using mosaik , allowing for zero mismatches between the reads and transcripts.
any parts of the transcripts for which mapped reads overlapped by less than 6 bp were considered hypothetical exon-exon junctions, and transcripts were split into “exons” at these positions.
using d. pseudoobscura gene models, it was confirmed that exon-exon junctions were correctly identified in a random test sample of transcripts.
next, both male and female genomic reads were mapped against the exons of the male muller c transcripts using mosaik, this time allowing for mismatches between the genomic reads and transcripts .
freebayes  was used to call snps in the mapping of male and female genomic reads versus muller c exons .
snps were called when  <dig> or more reads for the alternative alleles with a base quality of > <dig> were present at a given site, or if the coverage was  <dig> or less, the minimum number of reads supporting the alternative alleles was set to two.
based on the freebayes output, snps and indels detected by mapping genomic reads against the exons of transcripts were categorized into the following three classes :
<dig>  at a given position, male reads were polymorphic, and all female reads showed a variant that was different compared to the reference “exon” sequence .
these sites were assumed to already contain the neo-y variant in the transcript sequence.
<dig>  at a given position, both male and female reads showed a snp or indel; this could be caused by polymorphisms on the neo-x, sequencing error or by reads from paralogous locations being wrongly mapped to that particular position in both sexes; these positions were left unchanged .
<dig>  at a given position, male reads were polymorphic, whereas female reads showed no snp or indel .
such a site was considered to contain the neo-x variant in the transcript sequence and was modified to contain the neo-y variant, using a perl script.
genomic neo-y bac sequences  <cit>  were used to heuristically fine-tune the mosaik mapping parameters as well as the freebayes snp calling parameters, so as to minimize the neo-y versus bac sequence divergence.
our pipeline makes the implicit assumption that neo-x genomic reads can be mapped back to neo-y transcripts, which allows the identification of neo-x/neo-y diverged sites.
this assumption may be violated in regions where divergence is high enough for neo-x reads not to map, even if the homologous sequence is still present on the neo-x chromosome.
however, such regions are unlikely to affect the neo-y transcripts in the final assembly, since for regions of high divergence, trinity will assemble neo-x and neo-y transcripts separately.
our pipeline only modifies parts of a transcript where the neo-x and neo-y are similar enough for neo-x reads to align to the neo-y, and highly diverged regions will not be modified by our bioinformatics pipeline.
to test how this bioinformatic pipeline performs, we also mapped male and female reads against autosomal transcripts, and categorized variants in the same manner as for the neo-sex linked genes.
we find that  <dig> % of autosomal sites would fall into category  <dig>  and  <dig> % in category  <dig>  i.e.
the proportions of such sites relative to those for neo-sex-linked transcripts are 12% and 13%, respectively.
this suggests that about 12-13% of the sites on muller element c that appear as “male polymorphic and female homomorphic” may be caused by factors other than true neo-x/neo-y divergence.
we find that  <dig> % of autosomal sites would fall into category  <dig> .
since the proportions of category  <dig> snps are very similar between neo-sex chromosomes and autosomes, this is consistent with snps in this category being caused by mapping and sequencing errors, or remaining heterozygosity in the line, rather than true neo-x/neo-y divergence.
after modification of the male muller c transcripts, the male rna-seq reads were re-mapped against the modified transcripts, allowing for zero mismatches between reads and transcripts.
based on this mapping and the freebayes snp/indel output, transcript regions were kept for further analysis if they were supported by male rna-seq reads  and if they contained at least one variant distinguishing the neo-y from the neo-x.
requiring at least one neo-y-specific variant per contiguous transcript region has the consequence that genes that were deleted from the neo-y were removed from the neo-y assembly, as were neo-x/neo-y homomorphic genes and silenced neo-y genes.
using this approach alone, we missed male-specific neo-x transcripts because they were turned into neo-y transcripts; however, male-specific neo-x transcripts were assembled in an analogous way to the neo-y transcripts , except that this time positions were modified to carry the variant that was also found in females .
all transcripts were merged using cap <dig>  <cit> , and then scaffolded based on drosophila pseudoobscura known proteins using the stm method  <cit> .
all transcripts were assigned to flybase polypeptide ids, using a reciprocal best hit between the assembled transcripts and d. pseudoobscura flybase transcripts as a criterion for id assignment.
transcript level calculation
to calculate transcript abundance for the neo-x and neo-y separately, rna-seq reads from male carcass, whole body, testis and accessory gland were aligned to the d. miranda transcriptome assembly using mosaik and allowing for zero mismatches between the reads and the reference transcripts.
express  <cit>  probabilistically assigns reads to different alleles of the same gene – in this case, the neo-x and neo-y homologs – and was used to calculate transcript abundance levels  for the neo-sex chromosomes separately, as well as for genes on other muller elements.
sequence analysis of transcripts
neo-x and neo-y genes, as well as their autosomal d. pseudoobscura homologs were aligned using muscle  <cit> .
gblocks  <cit>  was used to remove poorly aligned segments of alignments, and the minimum length of a block was set to  <dig>  and the maximum number of consecutive non-conserved positions to four.
next, only alignments that contained all three sequences  were considered .
for these transcripts, distmat  was used to calculate the uncorrected divergence between the neo-x, neo-y and d. pseudoobscura.
to assess the functional properties of neo-y genes, the neo-y transcripts were aligned to the protein sequences of their d. pseudoobscura homologs using genewise  <cit> .
for the alignments, genewise allows for frame-shift mutations and pre-mature stop codons in the transcript sequences, and a perl program was used to extract the positions of these mutations within the neo-y transcripts.
for each neo-y transcript, the length of the 3′ utr was calculated by the distance between the first stop codon to the end of the transcript.
since neo-y expression is generally low, we also used the length of the 3′utr of the neo-x homolog as a proxy for utr length on the neo-y, assuming that the neo-y transcript might not be fully assembled if the read number supporting that part of the transcript was too low.
availability of supporting data
the data sets supporting the results of this article are available in the ncbi sequence read archive  repository, http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term=srr <dig>  http://www.ncbi.nlm.nih.gov/sra/?term= srr <dig>
RESULTS
the neo-sex transcriptome
the raw trinity assembly of the male transcriptome contained  <dig>  transcripts,  <dig>  of which showed sequence similarity  to muller element c of the d. miranda genomic assembly  <cit> .
it was manually confirmed that neo-sex transcripts in males were produced as chimeric transcripts by trinity, by comparing several transcripts to neo-y and neo-x bac clone sequences of ref.
<cit> .
the female trinity assembly contained  <dig>  transcripts,  <dig>  of which potentially neo-x linked.
<dig>  of the male muller c matching transcripts had exact matches with the neo-x assembly, suggesting that trinity assembled some neo-x transcripts correctly, in addition to chimeric and correctly assembled neo-y transcripts.
the  <dig>  potential muller c transcripts from males were divided into  <dig>  exonic sequences.
by mapping of genomic reads followed by snp calling, about  <dig> % of all sites in muller c transcripts from males were transformed into the corresponding neo-y variant.
all transcripts were scaffolded using the stm method and d. pseudoobscura protein sequences , which contains  <dig>  annotated polypeptides on muller c. the final d. miranda transcriptome assembly, which has been submitted to ncbi under the accession number galp <dig>  contained  <dig>  protein-coding genes, including  <dig>  from the neo-x,  <dig>  from the neo-y and  <dig>  from the autosomes/ancestral sex chromosomes.
thus, our de novo transcriptome captures a large fraction of genes present on the neo-sex chromosomes, with  <dig>  protein-coding genes being expressed from both the neo-x and neo-y,  <dig> transcripts only being detected from the neo-x, and  <dig> transcripts only from the neo-y.
n <dig> is a statistical measure of the average length of a set of sequences, defined as the length n for which 50% of all bases in the sequences are in a sequence of length l < n; n <dig> was  <dig>  bp for the neo-x;  <dig>  bp for the neo-y;  <dig>  bp for the autosomes and the older x chromosomes combined.
thus, neo-y transcripts were, on average, shorter compared to those from the neo-x or autosomes .
nevertheless, we obtain nearly as many transcripts from the neo-y as from the neo-x.
the median sequence divergence between the neo-y and neo-x was  <dig> % and between the neo-y and d. pseudoobscura  <dig> %, suggesting that the neo-sex chromosome system in d. miranda was formed shortly after the split of the two species, rather than at the time of speciation .
d. pseudoobscura
transcripts, respectively.
using genewise, neo-y transcripts were divided into those with putative functional orfs and those that are non-functional.
of the  <dig>  neo-y transcripts with no detectable frame-shift mutations or premature stop codons, the median length was  <dig>  – somewhat shorter than for transcripts annotated as containing disrupted orfs by frameshift mutations  or pre-mature stop codons  – combined,  <dig> putative non-functional neo-y genes, with a median length of  <dig>  bp .
genes containing nonsense or frameshift mutations might simply be longer, on average, because they present larger mutational targets, but notably, they are still transcribed throughout most of their length.
on the other hand, some knock-out mutations might have been missed if they were not part of the assembled transcript sequence.
indeed, compared to the genome assembly  <cit> , we detect fewer loss-of-function mutations among neo-y-linked genes ; however, we also find mutations in  <dig> genes that were annotated as “functional” in the genomic assembly .
partly, this may be due to the relatively high chance of missing a pseudogene in the genome assembly - which was estimated as 10%  <cit> .
in addition, we followed a different strategy of assigning flybase polypeptide ids to neo-y transcripts, and the two assemblies may differ with regards to which exact transcript is assigned to which fbpp id .
however, since all loss-of function mutations in the transcriptome assembly are supported by genomic as well as rna-seq reads, they are presumably real.
nonsense mediated decay
the median length of the 3′utr for the neo-x sequences was 346-bp, slightly shorter than the median utr length of d. pseudoobscura transcripts   <cit> .
in contrast, for neo-y genes without ptcs, the 3′utrs had a median length of only 112-bp, i.e.
neo-y 3'utrs were substantially shorter compared to the neo-x 3′utrs .
in part, this may be caused by trinity not being able to fully assemble neo-y utrs that are lowly expressed.
among the  <dig> neo-y genes with ptcs, median original utr length was only  <dig> -bp when counting from the ancestral stop codon to the end of the transcript, but it was increased to  <dig> -bp when counting from the first ptc to the end of the transcript.
accordingly, stop codons often occurred in the beginning of the transcript sequence, and appear distributed randomly across the coding sequence of a gene .
however, several lines of evidence suggest that nmd does not play a major role in down-regulating transcript levels of neo-y genes that contain ptcs.
first, fpkm-values were reduced for non-functional neo-y genes, irrespective of whether they contained a frame-shift mutation only without causing a ptc further down-stream of the transcript, or if they contained a ptc .
accordingly, most down-regulation of non-functional genes occurs irrespective of nmd, which suggest general mechanisms such as epigenetic modifications causing reduced expression of non-functional genes  <cit> , rather than down-regulation due to feedback in the cell after transcripts have been produced.
in line with this, neo-y genes that were not present in the de novo transcriptome assembly were often annotated as functional by ref.
<cit>  , and only 15% contained ptcs, i.e.
down-regulation by nmd cannot explain the lack of expression for the majority of these silenced transcripts.
furthermore, for transcripts containing ptcs, there was no detectable relationship between 3′utr length and transcript levels .
thus, our results suggest that nmd does not play a major role in down-regulating neo-y transcript levels in d. miranda if ptcs are present.
this is consistent with possibly abundant read-through of stop codons in drosophila, as indicated by protein-coding conservation 3′ of stop codons  <cit> .
in order to infer if translation proceeds past the pre-mature stop codons of neo-y-linked genes, direct measures of translation are necessary, such as ribosomal profiling.
interestingly, neo-x genes in d. miranda show significantly higher transcript levels if they contain utr lengths larger than the median .
if longer 3′utrs are indeed associated with higher levels of expression , this may counteract any effects of nmd.
gene loss and silencing on the neo-x chromosome
the availability of neo-y transcripts allowed us to identify genes that were lost from the neo-x since its split from the neo-y about  <dig>  my ago.
when mapping female genomic reads back to male transcripts, four muller c transcripts had zero female genomic read coverage, indicating that these genes are not present in the female genome.
there was also no polymorphism among male genomic reads mapped to these transcripts, further supporting the idea that only a single, neo-y-linked transcript was present in the genome.
a fifth gene  was identified as having an unusually high “divergence” from bac clone neo-y sequences of ref.
<cit> .
this gene turned out to be present in two copies in the neo-y genomic assembly and also in two copies on muller c in d. pseudoobscura, but only in one copy on the neo-x .
whole genome alignments of d. miranda, d. pseudoobscura and d. affinis <cit>  were used to verify the deletions on the neo-x, and to investigate the scale of the deletion.
in each case, neo-x deletions only affected single genes or parts of genes, suggesting that they were gene-specific and small in scale .
interestingly, the neo-y homologs of genes that were lost on the neo-x chromosome are highly expressed from the testis in d. miranda .
similarly, the  <dig> transcripts that were recovered from the neo-y and not the neo-x show testis-biased expression in d. miranda males .
expression levels of muller c genes in d. pseudoobscura can be used as a proxy to infer their ancestral function, before becoming sex-linked in d. miranda.
in the d. pseudoobscura transcriptome, homologs of only three of the five neo-x deleted genes could be retrieved; these, however, show almost no expression in females and moderate to high expression in male tissues , suggesting that they were ancestrally male-specific in function.
likewise, genes that have been silenced on the neo-x in the d. miranda lineage  are also male-biased in expression in d. pseudoobscura .
demasculinization  has been found on the ancestral x chromosome of drosophila, and on chromosome xr of the d. pseudoobscura lineage  <cit> .
our observations of both beginning down-regulation of male-specific transcripts as well as targeted gene loss of male-specific genes from the young neo-x of d. miranda demonstrate that these changes can occur rapidly within a short evolutionary time frame.
in contrast, the  <dig> genes that were recovered from the neo-x but not the neo-y  had generally low expression in d. miranda females and in both sexes of d. pseudoobscura , supporting the idea that lowly expressed genes are lost at an increased rate from the neo-y  <cit> .
mapping of male genomic reads against neo-x transcripts suggests that most silent genes on the neo-y are indeed still present in the genome, with only their transcription being suppressed, i.e.
their genomic coverage in males is only partially reduced .
indeed only  <dig> of the  <dig> genes with silent neo-y but expressed neo-x copies were annotated as deleted from the neo-y  <cit> .
CONCLUSIONS
evolved sex chromosomes have a highly diverged gene content  <cit> .
ancient y chromosomes have usually lost most of their ancestral genes, and the few remaining genes often have testis-specific functions  <cit> .
in contrast, old x chromosomes in drosophila have become demasculinized, containing a deficiency of genes with male-specific expression  <cit> .
the evolutionary events leading to such a difference in gene content, and their temporal dynamics, are little understood.
here, we show that both gene loss and gene silencing contribute to divergence in gene content between evolving sex chromosomes.
genes with male-specific function are preferentially silenced and lost on an evolving neo-x chromosome.
on the other hand, lowly expressed genes are lost quickly from the neo-y, while male-specific genes are more likely to be retained on the degenerating neo-y chromosome.
thus, masculinization of the y and demasculinization of the x can proceed very quickly after a new sex chromosome emerges, supporting the idea that sex-specific or sexually antagonistic selection plays a major role in the evolution of heteromorphic sex chromosomes  <cit> .
competing interests
the authors declare that they have no competing interests.
authors’ contributions
vbk built the bioinformatic pipeline, performed the analyses and wrote the paper; db conceived the study and contributed to the analysis and writing of the paper.
both authors read and approved the final manuscript.
supplementary material
additional file 1
venn diagram of the number of genes annotated as non-functional, i.e.
containing frame-shift mutations or ptcs, in the de novo transcriptome assembly and the genome assembly of <cit> .
this diagram is based on a total of  <dig>  neo-y genes that were allocated to the muller c element in both assemblies, expressed from the neo-y  and had flybase polypeptide ids assigned to them .
click here for file
additional file 2
detailed description of the neo-x deleted gene regions.
click here for file
additional file 3
genomic coverage of neo-x transcripts in males, for neo-x genes whose neo-y homologs are not transcribed and possibly deleted  or present in the neo-y transcriptome .
click here for file
acknowledgments
we thank chris ellison for the whole genome alignments of drosophila miranda, d. pseudoobscura and d. affinis.
the study was funded by grants from the national institutes of health  to db.

 the differentiation of an evolving pair of sex chromosomes involves many evolutionary steps, including changes in gene content on the x and y; gene loss and the fixation of deleterious mutations on the y; epigenetic modifications related to silencing of y-linked genes and the evolution of dosage compensation. the drosophila miranda neo-sex chromosome system is emerging as a useful resource to study the evolution of young sex chromosomes. d. miranda split from its sister species d. pseudoobscura only about my ago, and has sex chromosomes of different ages: muller a, the ancestral sex chromosome shared among all drosophila species including d. melanogaster, is at least my old; muller d, a neo-sex chromosome system of intermediate age and shared with d. pseudoobscura; muller c, specific to the d. miranda lineage. the genome sequence of d. miranda has recently been published, including a draft assembly of the neo-y chromosome. a major challenge in the assembly of this chromosome was the fact that the neo-y has already accumulated many repetitive sequences, which are difficult to assemble – similar to the ancestral y chromosome of drosophila which has not been assembled at all. in addition, the high level of sequence identity between the neo-x and neo-y especially in coding regions further complicates assemblies using next-generation sequencing data, since many short sequencing reads are identical from the neo-x and neo-y. to annotate the neo-y and study expression changes, diagnostic snps between the neo-x and neo-y were used to estimate the relative abundance of neo-sex linked transcripts. this, however, does not allow studying the structure of the neo-y transcripts per se. in particular, it is unknown how complete neo-y protein-coding genes are, especially with regards to the transcripts that contain premature stop codons or frameshift mutations. furthermore, the use of diagnostic snps for expression analysis will lead to an inherent mapping bias towards higher neo-x expression, since rna-seq reads can overlap with parts of transcripts that differ structurally between the neo-x and neo-y. using the neo-x sequence as a baseline for comparison, one can only study neo-y transcripts that are also present on the neo-x and show no major structural changes or indels. thus, there is considerable interest to study the neo-sex transcriptome directly, which we do here. in addition to providing insights into sex chromosome evolution, the neo-sex chromosomes of d. miranda also allow us to study other cellular phenomena, such as nonsense-mediated mrna decay. nmd is a cellular surveillance mechanism that degrades transcripts containing premature termination codons. destruction of transcripts containing ptcs by nmd prevents production of truncated, potentially harmful proteins that may interfere with normal cellular processes. the exact mechanism through which nmd recognized ptcs in drosophila is unclear, though the length of the 3′ utr has been implicated, whereas the presence of an exon-exon junction downstream of a stop codon is an initiator of nmd in mammals. given that a large number of neo-y transcripts contain ptcs, the d. miranda system makes an excellent case for testing if the length of the 3′utr leads to nmd in flies. a major challenge in constructing the neo-sex transcriptome is that coding regions on the neo-x and neo-y are very similar; based on the genome sequence, the estimated divergence at coding regions is only about %. accordingly, using illumina rna-seq reads for transcriptome assembly, not all reads can be assigned to the neo-x or neo-y unambiguously. this problem is for the most part irrelevant for the neo-x transcriptome, which can be assembled directly in females. in males, however, transcriptome assemblers fail to assemble the neo-x and neo-y into separate transcripts. instead, chimeric transcripts are produced that contain partial neo-x and partial neo-y sequence, connected by regions that are not differentiated between the neo-sex chromosomes. to obtain an assembly of the neo-y transcriptome, we developed a bioinformatic pipeline, making use of genomic read mapping against the neo-sex transcripts in males. methods assembly pipeline we used illumina genome analyzer ii paired-end 100-bp rna-seq reads of the inbred line msh to assemble the d. miranda transcriptome with trinity. first, low-quality sections of rna-seq reads and remaining adaptor sequence were discarded, and the minimum contig size of the trinity assembly was set to 200 bp. for the neo-y assembly, million rna-seq reads from four different samples were combined in a single assembly; similarly, to obtain neo-x transcripts, million rna-seq reads from virgin ovary, virgin female whole body and virgin female gonadectomized carcass were combined and assembled together. for the autosomes and old sex chromosomes, both male and female reads were combined for a single assembly. the d. pseudoobscura transcriptome was assembled de novo using combined virgin male and virgin female whole body samples. an overview of the pipeline to assemble the neo-y transcripts is given in figure  using blastn, trinity transcripts in males were compared to the d. miranda genome assembly. potential muller c transcripts were kept and modified to become neo-y transcripts, by introducing variants that distinguish the neo-x from the neo-y. for this means, male genomic reads were compared to the potential muller c transcripts and snp positions were identified. note that, in contrast to rna-seq reads, male genomic reads contain equal proportions of neo-x and neo-y variants in males. to avoid introducing erroneous variants at intron-exon junctions, male trinity transcripts were first divided up into exons. this was done by mapping male genomic reads against the male muller c transcripts using mosaik, allowing for zero mismatches between the reads and transcripts. any parts of the transcripts for which mapped reads overlapped by less than 6 bp were considered hypothetical exon-exon junctions, and transcripts were split into “exons” at these positions. using d. pseudoobscura gene models, it was confirmed that exon-exon junctions were correctly identified in a random test sample of transcripts. next, both male and female genomic reads were mapped against the exons of the male muller c transcripts using mosaik, this time allowing for mismatches between the genomic reads and transcripts. freebayes was used to call snps in the mapping of male and female genomic reads versus muller c exons. snps were called when or more reads for the alternative alleles with a base quality of > were present at a given site, or if the coverage was or less, the minimum number of reads supporting the alternative alleles was set to two. based on the freebayes output, snps and indels detected by mapping genomic reads against the exons of transcripts were categorized into the following three classes: at a given position, male reads were polymorphic, and all female reads showed a variant that was different compared to the reference “exon” sequence. these sites were assumed to already contain the neo-y variant in the transcript sequence. at a given position, both male and female reads showed a snp or indel; this could be caused by polymorphisms on the neo-x, sequencing error or by reads from paralogous locations being wrongly mapped to that particular position in both sexes; these positions were left unchanged. at a given position, male reads were polymorphic, whereas female reads showed no snp or indel. such a site was considered to contain the neo-x variant in the transcript sequence and was modified to contain the neo-y variant, using a perl script. genomic neo-y bac sequences were used to heuristically fine-tune the mosaik mapping parameters as well as the freebayes snp calling parameters, so as to minimize the neo-y versus bac sequence divergence. our pipeline makes the implicit assumption that neo-x genomic reads can be mapped back to neo-y transcripts, which allows the identification of neo-x/neo-y diverged sites. this assumption may be violated in regions where divergence is high enough for neo-x reads not to map, even if the homologous sequence is still present on the neo-x chromosome. however, such regions are unlikely to affect the neo-y transcripts in the final assembly, since for regions of high divergence, trinity will assemble neo-x and neo-y transcripts separately. our pipeline only modifies parts of a transcript where the neo-x and neo-y are similar enough for neo-x reads to align to the neo-y, and highly diverged regions will not be modified by our bioinformatics pipeline. to test how this bioinformatic pipeline performs, we also mapped male and female reads against autosomal transcripts, and categorized variants in the same manner as for the neo-sex linked genes. we find that % of autosomal sites would fall into category and % in category i.e. the proportions of such sites relative to those for neo-sex-linked transcripts are 12% and 13%, respectively. this suggests that about 12-13% of the sites on muller element c that appear as “male polymorphic and female homomorphic” may be caused by factors other than true neo-x/neo-y divergence. we find that % of autosomal sites would fall into category. since the proportions of category snps are very similar between neo-sex chromosomes and autosomes, this is consistent with snps in this category being caused by mapping and sequencing errors, or remaining heterozygosity in the line, rather than true neo-x/neo-y divergence. after modification of the male muller c transcripts, the male rna-seq reads were re-mapped against the modified transcripts, allowing for zero mismatches between reads and transcripts. based on this mapping and the freebayes snp/indel output, transcript regions were kept for further analysis if they were supported by male rna-seq reads and if they contained at least one variant distinguishing the neo-y from the neo-x. requiring at least one neo-y-specific variant per contiguous transcript region has the consequence that genes that were deleted from the neo-y were removed from the neo-y assembly, as were neo-x/neo-y homomorphic genes and silenced neo-y genes. using this approach alone, we missed male-specific neo-x transcripts because they were turned into neo-y transcripts; however, male-specific neo-x transcripts were assembled in an analogous way to the neo-y transcripts, except that this time positions were modified to carry the variant that was also found in females. all transcripts were merged using cap, and then scaffolded based on drosophila pseudoobscura known proteins using the stm method. all transcripts were assigned to flybase polypeptide ids, using a reciprocal best hit between the assembled transcripts and d. pseudoobscura flybase transcripts as a criterion for id assignment. transcript level calculation to calculate transcript abundance for the neo-x and neo-y separately, rna-seq reads from male carcass, whole body, testis and accessory gland were aligned to the d. miranda transcriptome assembly using mosaik and allowing for zero mismatches between the reads and the reference transcripts. express probabilistically assigns reads to different alleles of the same gene – in this case, the neo-x and neo-y homologs – and was used to calculate transcript abundance levels for the neo-sex chromosomes separately, as well as for genes on other muller elements. sequence analysis of transcripts neo-x and neo-y genes, as well as their autosomal d. pseudoobscura homologs were aligned using muscle. gblocks was used to remove poorly aligned segments of alignments, and the minimum length of a block was set to and the maximum number of consecutive non-conserved positions to four. next, only alignments that contained all three sequences were considered. for these transcripts, distmat was used to calculate the uncorrected divergence between the neo-x, neo-y and d. pseudoobscura. to assess the functional properties of neo-y genes, the neo-y transcripts were aligned to the protein sequences of their d. pseudoobscura homologs using genewise. for the alignments, genewise allows for frame-shift mutations and pre-mature stop codons in the transcript sequences, and a perl program was used to extract the positions of these mutations within the neo-y transcripts. for each neo-y transcript, the length of the 3′ utr was calculated by the distance between the first stop codon to the end of the transcript. since neo-y expression is generally low, we also used the length of the 3′utr of the neo-x homolog as a proxy for utr length on the neo-y, assuming that the neo-y transcript might not be fully assembled if the read number supporting that part of the transcript was too low. availability of supporting data the data sets supporting the results of this article are available in the ncbi sequence read archive repository, http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term=srr http://www.ncbi.nlm.nih.gov/sra/?term= srr  the neo-sex transcriptome the raw trinity assembly of the male transcriptome contained transcripts, of which showed sequence similarity to muller element c of the d. miranda genomic assembly. it was manually confirmed that neo-sex transcripts in males were produced as chimeric transcripts by trinity, by comparing several transcripts to neo-y and neo-x bac clone sequences of ref.. the female trinity assembly contained transcripts, of which potentially neo-x linked. of the male muller c matching transcripts had exact matches with the neo-x assembly, suggesting that trinity assembled some neo-x transcripts correctly, in addition to chimeric and correctly assembled neo-y transcripts. the potential muller c transcripts from males were divided into exonic sequences. by mapping of genomic reads followed by snp calling, about % of all sites in muller c transcripts from males were transformed into the corresponding neo-y variant. all transcripts were scaffolded using the stm method and d. pseudoobscura protein sequences, which contains annotated polypeptides on muller c. the final d. miranda transcriptome assembly, which has been submitted to ncbi under the accession number galp contained protein-coding genes, including from the neo-x, from the neo-y and from the autosomes/ancestral sex chromosomes. thus, our de novo transcriptome captures a large fraction of genes present on the neo-sex chromosomes, with protein-coding genes being expressed from both the neo-x and neo-y, transcripts only being detected from the neo-x, and transcripts only from the neo-y. n is a statistical measure of the average length of a set of sequences, defined as the length n for which 50% of all bases in the sequences are in a sequence of length l < n; n was  bp for the neo-x;  bp for the neo-y;  bp for the autosomes and the older x chromosomes combined. thus, neo-y transcripts were, on average, shorter compared to those from the neo-x or autosomes. nevertheless, we obtain nearly as many transcripts from the neo-y as from the neo-x. the median sequence divergence between the neo-y and neo-x was % and between the neo-y and d. pseudoobscura %, suggesting that the neo-sex chromosome system in d. miranda was formed shortly after the split of the two species, rather than at the time of speciation. d. pseudoobscura transcripts, respectively. using genewise, neo-y transcripts were divided into those with putative functional orfs and those that are non-functional. of the neo-y transcripts with no detectable frame-shift mutations or premature stop codons, the median length was – somewhat shorter than for transcripts annotated as containing disrupted orfs by frameshift mutations or pre-mature stop codons – combined, putative non-functional neo-y genes, with a median length of  bp. genes containing nonsense or frameshift mutations might simply be longer, on average, because they present larger mutational targets, but notably, they are still transcribed throughout most of their length. on the other hand, some knock-out mutations might have been missed if they were not part of the assembled transcript sequence. indeed, compared to the genome assembly, we detect fewer loss-of-function mutations among neo-y-linked genes; however, we also find mutations in genes that were annotated as “functional” in the genomic assembly. partly, this may be due to the relatively high chance of missing a pseudogene in the genome assembly - which was estimated as 10%. in addition, we followed a different strategy of assigning flybase polypeptide ids to neo-y transcripts, and the two assemblies may differ with regards to which exact transcript is assigned to which fbpp id. however, since all loss-of function mutations in the transcriptome assembly are supported by genomic as well as rna-seq reads, they are presumably real. nonsense mediated decay the median length of the 3′utr for the neo-x sequences was 346-bp, slightly shorter than the median utr length of d. pseudoobscura transcripts. in contrast, for neo-y genes without ptcs, the 3′utrs had a median length of only 112-bp, i.e. neo-y 3'utrs were substantially shorter compared to the neo-x 3′utrs. in part, this may be caused by trinity not being able to fully assemble neo-y utrs that are lowly expressed. among the neo-y genes with ptcs, median original utr length was only -bp when counting from the ancestral stop codon to the end of the transcript, but it was increased to -bp when counting from the first ptc to the end of the transcript. accordingly, stop codons often occurred in the beginning of the transcript sequence, and appear distributed randomly across the coding sequence of a gene. however, several lines of evidence suggest that nmd does not play a major role in down-regulating transcript levels of neo-y genes that contain ptcs. first, fpkm-values were reduced for non-functional neo-y genes, irrespective of whether they contained a frame-shift mutation only without causing a ptc further down-stream of the transcript, or if they contained a ptc. accordingly, most down-regulation of non-functional genes occurs irrespective of nmd, which suggest general mechanisms such as epigenetic modifications causing reduced expression of non-functional genes, rather than down-regulation due to feedback in the cell after transcripts have been produced. in line with this, neo-y genes that were not present in the de novo transcriptome assembly were often annotated as functional by ref., and only 15% contained ptcs, i.e. down-regulation by nmd cannot explain the lack of expression for the majority of these silenced transcripts. furthermore, for transcripts containing ptcs, there was no detectable relationship between 3′utr length and transcript levels. thus, our results suggest that nmd does not play a major role in down-regulating neo-y transcript levels in d. miranda if ptcs are present. this is consistent with possibly abundant read-through of stop codons in drosophila, as indicated by protein-coding conservation 3′ of stop codons. in order to infer if translation proceeds past the pre-mature stop codons of neo-y-linked genes, direct measures of translation are necessary, such as ribosomal profiling. interestingly, neo-x genes in d. miranda show significantly higher transcript levels if they contain utr lengths larger than the median. if longer 3′utrs are indeed associated with higher levels of expression, this may counteract any effects of nmd. gene loss and silencing on the neo-x chromosome the availability of neo-y transcripts allowed us to identify genes that were lost from the neo-x since its split from the neo-y about my ago. when mapping female genomic reads back to male transcripts, four muller c transcripts had zero female genomic read coverage, indicating that these genes are not present in the female genome. there was also no polymorphism among male genomic reads mapped to these transcripts, further supporting the idea that only a single, neo-y-linked transcript was present in the genome. a fifth gene was identified as having an unusually high “divergence” from bac clone neo-y sequences of ref.. this gene turned out to be present in two copies in the neo-y genomic assembly and also in two copies on muller c in d. pseudoobscura, but only in one copy on the neo-x. whole genome alignments of d. miranda, d. pseudoobscura and d. affinis were used to verify the deletions on the neo-x, and to investigate the scale of the deletion. in each case, neo-x deletions only affected single genes or parts of genes, suggesting that they were gene-specific and small in scale. interestingly, the neo-y homologs of genes that were lost on the neo-x chromosome are highly expressed from the testis in d. miranda. similarly, the transcripts that were recovered from the neo-y and not the neo-x show testis-biased expression in d. miranda males. expression levels of muller c genes in d. pseudoobscura can be used as a proxy to infer their ancestral function, before becoming sex-linked in d. miranda. in the d. pseudoobscura transcriptome, homologs of only three of the five neo-x deleted genes could be retrieved; these, however, show almost no expression in females and moderate to high expression in male tissues, suggesting that they were ancestrally male-specific in function. likewise, genes that have been silenced on the neo-x in the d. miranda lineage are also male-biased in expression in d. pseudoobscura. demasculinization has been found on the ancestral x chromosome of drosophila, and on chromosome xr of the d. pseudoobscura lineage. our observations of both beginning down-regulation of male-specific transcripts as well as targeted gene loss of male-specific genes from the young neo-x of d. miranda demonstrate that these changes can occur rapidly within a short evolutionary time frame. in contrast, the genes that were recovered from the neo-x but not the neo-y had generally low expression in d. miranda females and in both sexes of d. pseudoobscura, supporting the idea that lowly expressed genes are lost at an increased rate from the neo-y. mapping of male genomic reads against neo-x transcripts suggests that most silent genes on the neo-y are indeed still present in the genome, with only their transcription being suppressed, i.e. their genomic coverage in males is only partially reduced. indeed only of the genes with silent neo-y but expressed neo-x copies were annotated as deleted from the neo-y.  evolved sex chromosomes have a highly diverged gene content. ancient y chromosomes have usually lost most of their ancestral genes, and the few remaining genes often have testis-specific functions. in contrast, old x chromosomes in drosophila have become demasculinized, containing a deficiency of genes with male-specific expression. the evolutionary events leading to such a difference in gene content, and their temporal dynamics, are little understood. here, we show that both gene loss and gene silencing contribute to divergence in gene content between evolving sex chromosomes. genes with male-specific function are preferentially silenced and lost on an evolving neo-x chromosome. on the other hand, lowly expressed genes are lost quickly from the neo-y, while male-specific genes are more likely to be retained on the degenerating neo-y chromosome. thus, masculinization of the y and demasculinization of the x can proceed very quickly after a new sex chromosome emerges, supporting the idea that sex-specific or sexually antagonistic selection plays a major role in the evolution of heteromorphic sex chromosomes. competing interests the authors declare that they have no competing interests. authors’ contributions vbk built the bioinformatic pipeline, performed the analyses and wrote the paper; db conceived the study and contributed to the analysis and writing of the paper. both authors read and approved the final manuscript. supplementary material additional file 1 venn diagram of the number of genes annotated as non-functional, i.e. containing frame-shift mutations or ptcs, in the de novo transcriptome assembly and the genome assembly of. this diagram is based on a total of neo-y genes that were allocated to the muller c element in both assemblies, expressed from the neo-y and had flybase polypeptide ids assigned to them. click here for file additional file 2 detailed description of the neo-x deleted gene regions. click here for file additional file 3 genomic coverage of neo-x transcripts in males, for neo-x genes whose neo-y homologs are not transcribed and possibly deleted or present in the neo-y transcriptome. click here for file acknowledgments we thank chris ellison for the whole genome alignments of drosophila miranda, d. pseudoobscura and d. affinis. the study was funded by grants from the national institutes of health to db.
thus, evolving sex chromosomes can rapidly shift their gene content or patterns of gene expression in response to their sex-biased transmission, supporting the idea that sex-specific or sexually antagonistic selection plays a major role in the evolution of heteromorphic sex chromosomes. we find evidence of gene loss of male-specific genes on the neo-x chromosome, and transcriptional silencing of testis-specific genes from the neo-x. the drosophila miranda neo-sex chromosome system is a useful resource for studying recently evolved sex chromosomes. here, we de novo assemble the transcriptome of d. miranda using rna-seq data from several male and female tissues, and develop a bioinformatic pipeline to separately reconstruct neo-x and neo-y transcripts. we find that 24% of expressed neo-y transcripts harbor nonsense mutation within their open reading frames, yet most non-functional neo-y genes are expressed throughout all of their length. nonsense mediated decay has been implicated to degrade transcripts containing pre-mature termination codons in drosophila, but rampant description of neo-y genes with pre-mature stop codons suggests that it does not play a major role in down-regulating transcripts from the neo-y. neo-y transcripts are generally shorter than their homologous neo-x transcripts and expressed at lower levels. however, the neo-y genomic assembly is fragmented due to the accumulation of repetitive sequence.
 the drosophila miranda neo-sex chromosome system is a useful resource for studying recently evolved sex chromosomes. however, the neo-y genomic assembly is fragmented due to the accumulation of repetitive sequence. furthermore, the separate assembly of the neo-x and neo-y chromosomes into genomic scaffolds has proven to be difficult, due to their low level of sequence divergence, which in coding regions is about %. here, we de novo assemble the transcriptome of d. miranda using rna-seq data from several male and female tissues, and develop a bioinformatic pipeline to separately reconstruct neo-x and neo-y transcripts.  we obtain transcripts from the neo-x and from the neo-y. neo-y transcripts are generally shorter than their homologous neo-x transcripts and expressed at lower levels. we find that 24% of expressed neo-y transcripts harbor nonsense mutation within their open reading frames, yet most non-functional neo-y genes are expressed throughout all of their length. we find evidence of gene loss of male-specific genes on the neo-x chromosome, and transcriptional silencing of testis-specific genes from the neo-x.  nonsense mediated decay has been implicated to degrade transcripts containing pre-mature termination codons in drosophila, but rampant description of neo-y genes with pre-mature stop codons suggests that it does not play a major role in down-regulating transcripts from the neo-y. loss or transcriptional down-regulation of genes from the neo-x with male-biased function provides evidence for beginning demasculinization of the neo-x. thus, evolving sex chromosomes can rapidly shift their gene content or patterns of gene expression in response to their sex-biased transmission, supporting the idea that sex-specific or sexually antagonistic selection plays a major role in the evolution of heteromorphic sex chromosomes. sex chromosomesdrosophilatranscriptomegene lossnonsense mediated decay 