Transcriptome assemblies for studying sex-biased gene expression in the guppy, Poecilia reticulata
© Sharma et al.; licensee BioMed Central Ltd. 2014
Received: 6 January 2014
Accepted: 9 May 2014
Published: 26 May 2014
Skip to main content
© Sharma et al.; licensee BioMed Central Ltd. 2014
Received: 6 January 2014
Accepted: 9 May 2014
Published: 26 May 2014
Sexually dimorphic phenotypes are generally associated with differential gene expression between the sexes. The study of molecular evolution and genomic location of these differentially expressed, or sex-biased, genes is important for understanding inter-sexual divergence under sex-specific selection pressures. Teleost fish provide a unique opportunity to examine this divergence in the presence of variable sex-determination mechanisms of recent origin. The guppy, Poecilia reticulata, displays sexual dimorphism in size, ornaments, and behavior, traits shaped by natural and sexual selection in the wild.
To gain insight into molecular mechanisms underlying the guppy’s sexual dimorphism, we assembled a reference transcriptome combining genome-independent as well as genome-guided assemblies and analyzed sex-biased gene expression between different tissues of adult male and female guppies. We found tissue-associated sex-biased expression of genes related to pigmentation, signal transduction, and spermatogenesis in males; and growth, cell-division, extra-cellular matrix organization, nutrient transport, and folliculogenesis in females. While most sex-biased genes were randomly distributed across linkage groups, we observed accumulation of ovary-biased genes on the sex linkage group, LG12. Both testis-biased and ovary-biased genes showed a significantly higher rate of non-synonymous to synonymous substitutions (d N /d S ) compared to unbiased genes. However, in somatic tissues only female-biased genes, including those co-expressed in multiple tissues, showed elevated ratios of non-synonymous substitutions.
Our work identifies a set of annotated gene products that are candidate factors affecting sexual dimorphism in guppies. The differential genomic distribution of gonad-biased genes provides evidence for sex-specific selection pressures acting on the nascent sex chromosomes of the guppy. The elevated rates of evolution of testis-biased and female-biased genes indicate differing evolution under distinct selection pressures on the reproductive versus non-reproductive tissues.
In sexually reproducing species, males and females evolve differently due to differing regimes of natural and sexual selection [1–3]. Nonetheless, the evolution of sexually dimorphic traits within a species is constrained because most of the genome is shared between males and females. Therefore, the development of sex-specific traits is thought to be predominantly accomplished by sex-specific gene expression [4–8]. Quantitative analyses of complementary DNA (cDNA) from male and female tissues of mice (Mus musculus) , zebrafish (Danio rerio) , birds (chicken (Gallus gallus)  and turkey (Meleagris gallopavo) ), and insects (Drosophila species and Bombyx mori[13, 14]) have shown that a significant fraction of autosomal genes are differentially expressed between the sexes in their reproductive as well as non-reproductive tissues. This suggests that sex-biased gene expression contributes to sexually dimorphic phenotypes and sex-biased genes may evolve differently under selection pressures acting on the sexual phenotypes . Research on evolutionary properties of sex-biased genes has shown accelerated rates of coding sequence changes in reproduction-related male-biased genes [16, 17]. This is primarily attributed to greater sexual selection on males than females. Elevated nucleotide substitution rates of sex-biased genes expressed in somatic and reproductive tissues may also occur due to relaxed selection on non-pleiotropic tissue-specific genes [18, 19]. Sex-biased genes also show non-random genomic distribution with X- or Z-linkage [20, 21] that can arise due to differential selection on the hemizygous sex chromosome [6, 7].
So far, sex-biased gene expression has mainly been explored in species with well-differentiated sex chromosomes, while species with young or undifferentiated sex-chromosome systems are just beginning to be studied [22–24]. In this regard, teleost fish with their spectacular diversity of sex determination mechanisms and a large repertoire of duplicated genes provide largely unexplored resources to study sexual dimorphism resulting from sex-biased and sex-limited gene expression . Among teleosts, members of the family Poeciliidae are known to have multiple sex determination systems [26, 27] and are characterized by highly variable sexually dimorphic traits including color patterns, body size, genital morphology, and mating behavior [28–30]. The guppy (Poecilia reticulata) was one of the first vertebrates where XY sex determination and Y-linked inheritance of sexually selected traits were demonstrated . Sexual dimorphism in guppies is characterized by male-specific color patterns that attract females but are disadvantageous in the presence of predators [32–35]. These male-advantageous traits are believed to have co-evolved with female mate-choice preferences [36, 37]. The guppy also displays sexual size dimorphism. Female guppies grow throughout their life, whereas males slow down their growth during maturation . Male and female guppies also display behavioral differences in the amount of time spent mating, foraging, shoaling, and avoiding predators [32, 39–44].
While the evolutionary ecology of the guppy’s sexual dimorphism has been well studied with respect to heredity and adaptation, the molecular mechanisms governing this dimorphism are largely unidentified. Recently, using a high-density linkage map, quantitative trait loci (QTL) influencing male size, shape, and color traits were mapped to several sex-linked and autosomal loci . Nevertheless, in order to understand the contribution of sex-biased gene regulation to sexually dimorphic phenotypes, a genome-wide comparison of gene expression in sexually dimorphic tissues is required.
Current transcriptome resources of the guppy include a database of Sanger-sequenced expressed sequence tags (EST) and a more recent 454 sequenced transcriptome, that together correspond to roughly 9,000 unique transcripts from embryos and adult guppies originating from several different populations [46, 47].
Here we extend these resources by assembling a reference transcriptome using high depth Illumina sequencing. We used cDNA from multiple tissues from embryos and adults from a single guppy population, thereby minimizing population-specific effects in phenotypic variations and X- and Y-linkage [48, 49]. We then combined the predicted coding sequences from both genome-independent and genome-guided assemblers. The merged reference comprises expressed sequences from embryos and differentiated adult tissues. The transcriptome reference was then used to identify genes with either male- or female-biased expression in three tissues with phenotypic sexual dimorphism in the adult guppy. These included two somatic tissues (brain and tail) and the gonads. Furthermore, by examining sex-biased genes we determined whether i) the expression bias in adult guppy tissues reflects the morphological and physiological differences between the sexes; ii) there is non-random genomic distribution of these genes; and iii) they show signatures of relaxed selection when compared to unbiased genes, as hypothesized for genes subject to sexual selection.
Comparison of transcriptome assembled with genome-guided and genome-independent assemblers
Trinity: Genome-independent assembly (GIA)
Cufflinks: Genome-guided assembly (GGA)
Total length (bp)
Length with longest isoforms per locus (bp)
No. of transfrags
No. of transcripts (Unique components/gene groups)
Mean length (bp)
Longest contig (bp)
Overall mapping (%)
Concordant and unique mapping (%)
Total no. of ORFs
No. of complete ORFs
Mean length ORF (bp)
Longest ORF (bp)
Total length of assembly with CDSs only (bp)
Against guppy (GGA against GIA or vice-versa)
Orthologs in only one assembly
The genome-guided and genome-independent assemblies were compared using read-based, length-based, and annotation-based metrics. We compared the i) percent of reads remapped to the transcriptome (completeness); ii) the percent of correctly oriented mapped read pairs (accuracy); iii) total length of assembly and mean length of assembled transcripts (contiguity); iv) number and length of predicted ORFs, and v) number of orthologs identified using reciprocal Blastp against other validated protein sequence databases (Table 1).
Mapping the RNA-seq reads to each assembly we found that the genome-independent assembly incorporated a larger number of correctly oriented read pairs as compared to the genome-guided assembly (Table 1). On the other hand, the genome-guided assembly was more contiguous with longer transcripts, greater number of ORFs, and substantially more full-length ORFs (Table 1).
By examining the number of single-copy orthologs identified from comparing translated coding sequences (CDS) of the guppy against other teleost, human, and mouse protein sequence databases, we identified a greater number of orthologs in the genome-guided assembly than in the genome-independent assembly (Figure 1A, B, Table 1). The total number of orthologs found between guppy and other species was related to divergence between the species (with the exception of medaka, Oryzias latipes, possibly due to the smaller size of the medaka protein database) (Figure 1B). We identified 24,020 reciprocal best-blast hits shared between the genome-guided and genome-independent assemblies (Figure 1A, B, and Table 1). For approximately half (12,006) of these overlapping sequences, orthologous protein sequences were identified in other vertebrates. An additional 11,721 vertebrate orthologs were identified from only one of the two assemblies (Figure 1A, B, and Table 1). In addition to the identified orthologs, 30-40% of the remaining translated CDS aligned (alignment length > 50 amino acids) with significant sequence similarity (E-value < 1 × 10−20) to protein coding sequences of the other vertebrates (Table 1).
We merged the CDS predicted from transcripts of genome-guided and genome-independent assemblies to obtain a single comprehensive reference combining advantages from both assembly methods. This final dataset consisted of 74,567 sequences, hereafter referred to as the guppy reference transcriptome (GRT) (Figure 1C). In total, 30,643 (41.1% of the GRT) sequences showed significant alignment (Blastx E-value < 1 × 10−15) to 22,780 unique protein sequences in the NCBI non-redundant protein database (NR) . Out of these, 17,931 were annotated with functional categories (Gene Ontology: GO terms) (Figure 1C). A complete list of the best-blast hits and GO annotations is given in Additional file 2: Table S2 and Additional file 3: Table S3. A total of 73,518 sequences could be aligned to assembled scaffolds from the female genome. Of these, 67,882 aligned to scaffolds that were assigned to guppy linkage groups. All the sequences that did not align (1044) were from the genome-independent assembly and of these 693 (66.4%) could be aligned to the genome of the closely related platyfish, Xiphophorus maculatus (data not shown).
Differentially expressed genes between males and females in brain, gonad, and tail tissue
No. of expressed genes
Sex-biased genes (%)
Male-biased genes (%)
Female-biased genes (%)
> 1.2 fold
> Median fold
> 1.2 fold
> Median fold
> 1.2 fold
> Median fold
> 1.2 fold
> Median fold
> 1.2 fold
> Median fold
> 1.2 fold
> Median fold
Brain, Tail, Gonad
> 1.2 fold
> Median fold
Enriched Gene Ontology (GO) terms in male-bias and female-bias genes
(A) Brain tissue
Male: differentially expressed: 702; best Blastx hits in NR database: 420; genes with GO terms: 237
Ionotropic glutamate receptor signaling (GO:0035235: BP)
GRIK5 (2of2); GRIK4; GRIN2A (1 of 2); Glutamate [NMDA] receptor subunit zeta-1
Ion transmembrane transport (GO:0034220: BP)
Potassium channels: KCNJ3; KCNJ11 (1of2); KCNH1; KCNA1 (2 of2), Calcium channels: CACNA1I (4of4); Ryanodine receptor (RYR2 (1 of 3)), Sodium channel: SCN8A(2of2)
Regulation of cell development (GO:0060284: BP)
Lens calpain-3 (CAPN3); distal-less homeobox gene-1a (DLX1a); protein tyrosine phosphatase receptor-d (PTPRD); retinal cadherin-4 (CDH4 (1 of 2))
Cerebellar Purkinje cell differentiation (GO:0021702: BP)
LIM/homeobox protein (LHX1); voltage-dependent P/Q-type calcium channel subunit alpha-1A (CACNA1A) ; LIM domain binding-1 (LDB1(2 of 2))
Integral to membrane (GO:0016021: CC)
Hypocretin (orexin) receptor-1 (HCRTR1); multidrug resistance-associated protein 9; coiled-coil domain containing 149 CCDC149 (2 of 2); protocadherin 15b; gamma-aminobutyric acid A receptor, alpha-3 (GABRA3)
Female: differentially expressed: 1764; best Blastx hits in NR database: 1596; genes with GO terms: 955
Endothelial cell migration (GO:0043542: BP)
Sushi-repeat containing protein, X-linked2 (SRPX2); angiopoietin1 (ANGPT1); myosin heavy chain 9,non-muscle (MYH9(2 of 2))
Gonad development (GO:0008406: BP)
WNT10A; WNT5A; secreted frizzled related proteins (SFRP1, SFRP5) ; TGFB2; Phospholipase A2 groupIVa (PLAG4A);
Immune response (GO:0006955: BP)
WNT5A; TGFB2; kelch-like protein 6 (KLHL6); MHC classI-E (HLA-E); complement component (C6); MHC class II invariant chain (CD74(1of2)); chemokine(C-X-C motif) ligand12 (CXCL12(2of2))
Integrin-mediated signaling pathway (GO:0007229: BP)
Integrin b1 binding protein1 (ITGB1BP1); Integrin, alpha10 (ITGA10);Integrin beta (ITGBL1, ITGB3a); nicotinamide riboside kinase-2 (NMRK2(2of2))
DNA-dependent DNA replication initiation (GO:0006270: BP)
Minichromosome maintenance complex components (MCM2,MCM3,MCM4,MCM5,MCM6)
proteinaceous extracellular matrix (GO:0005578: CC)
ADAM metallopeptidase with thrombospondin family members (ADAMTS12, ADAMTS15); NID1 (1of2); FBN3; MFAP2; CYR61 (2of2); COL4A1; COL10A1; COL11A1; WNT11; matrix metalloproteinases (MMP12, MMP14,MMP2 )
(B) Tail tissue
Male: differentially expressed: 755; best Blastx hits in NR database: 635; genes with GO terms: 404
Neuropeptide signaling pathway (GO:0007218: BP)
Tachykinin, precursor 1 (TAC1); prepronociceptin PNOC (1 of 2); secretogranin V (7B2like) (SCG5); brain-specific angiogenesis inhibitor-3 (BAI3)
Neurotransmitter transport (GO:0006836: BP)
Syntaxin1B (STX1B); solute carrier family (SLC6A2, SLC6A5); syntaxin binding protein 1 STXBP1 (1 of 2)
Locomotory behavior (GO:0007626: BP)
Glycine receptor subunit beta (GLR-b2); choline O-acetyltransferase CHAT (2 of 2); astrotactin 1 (ASTN1)
5,6-dihydroxyindole-2-carboxylic acid oxidase TYRP1 (1 of 2, 2 of 2); Tyrosinase TYR (2 of 2); synaptotagmin-like-2 (SYTL2); adenosine deaminase CECR1
*Pigment biosynthetic process (GO:0046148: BP)
Insulin receptor binding (GO:0005158: MF)
Sorbin and SH3 domain containing 1 SORBS1; DOK7 (1 of 2); growth factor receptor-bound protein 10 (GRB10)
Female: differentially expressed: 705; best Blastx hits in NR database: 616; genes with GO terms: 387
Glycolysis (GO:0006096: BP)
lactate dehydrogenase A (LDHA); phosphoglycerate kinase 1 (PGK1); pyruvate kinase, muscle (PKM (1 of 2)); glucose-6-phosphate isomerase (GPI (2 of 2)); 2,3-bisphosphoglycerate mutase (BPGM)
DNA-dependent DNA replication (GO:0006261: BP)
Topoisomerase2a (TOP2a); primase, DNA, polypeptide 1 (PRIM1); tonsoku-like, DNA repair protein (TONSL); minichromosome maintenance complex components (MCM2, MCM4, MCM5, MCM6); polymerase (DNA directed), epsilon2 (POLE2)
Mitosis (GO:0007067: BP)
Aurora kinase C AURKC; cyclin-dependent kinase 1 CDK1 ; spindle apparatus coiled-coil protein 1 SPDL1; cyclin B1 CCNB1; non-SMC condensin I complex, subunitH, subunitD2 (NCAPH, NCAPD2); checkpoint kinase-1 (CHEK1)
Proteinaceous extracellular matrix (GO:0005578: CC)
Secreted protein, acidic, cysteine-rich (osteonectin) (SPARC); versican VCAN VCAN (2 of 2); tenascin(TNC (1 of 2)); collagens (COL11A1a, COL11A1b, COL27A1 (2 of 2))
(C) Gonad Tissue
Male: Differentially expressed: 4891; Best Blastx hits in NR database: 3879; Genes with GO terms: 2033
Cilium assembly (GO:0042384: BP)
Radial spoke head 9 homolog (RSPH9); ARPC4-TTLL3 readthrough; forkhead box J1 (FOXJ1); Transmembrane proteins (TMEM237, TMEM17, TMEM231); B9 protein domain 1 (B9D1); coiled-coil and C2 domain containing-2A (CC2D2A)
Spermatogenesis (GO:0007283: BP)
Kelch-like family member 10 KLHL10 (3 of 3); rhophilin associated tail protein 1-like (ROPN1L); phosphate cytidylyltransferase 1, choline-b ( PCYT1B (1 of 2)); Outer dense fiber protein 2/Cenexin (ODF2)
Meiosis I (GO:0007127: BP)
MutS homolog 5 (E. coli) MSH5; DMC1 dosage suppressor of mck1 homolog, meiosis-specific homologous recombination (DMC1); HORMA domain containing 1 HORMAD1; cyclin B1 interacting protein 1, E3 ubiquitin protein ligase CCNB1IP1
Fertilization (GO:0009566: BP)
KLHL10 (3 of 3); spectrin, beta, non-erythrocytic 4 SPTBN4 (1 of 2); glycine receptor, beta GLRB (1 of 2); polypyrimidine tract-binding protein-1 (PTPB1)
Female: Differentially expressed: 5163; Best Blastx hits in NR database: 4577; Genes with GO terms: 2847
Regulation of BMP signaling pathway (GO:0030510: BP)
Forkhead box H1 (FOXH1); Noggin (NOG); WNT5a; activin A receptor, type I (ACVR1); follistatin-like 1 (FSTL1); bone morphogenetic protein 15 (BMP15); growth differentiation factor 9 (Gdf9)
Fibroblast growth factor receptor signaling (GO:0008543: BP)
Kinesin family member-16Bb (KIF16bb); sal-like-4 (SALL4); serine threonine-protein phosphatase 2a (PP2a); sal-like 1 (SALL1); catenin, beta 1 (CTNNB1); fibroblast growth factor receptor-1(FGFR1 (2 of 2)); sprouty homolog 1, antagonist of FGF signaling (SPRY1); FGF20, FGF13, FGF 10
Focal adhesion (GO:0005925: CC)
Filamin A alpha FLNA (2 of 2); ezrin EZR (1 of 2); Rho guanine nucleotide exchange factor-7 (ARHGEF7 (al 3 paralogs)); syndecan-4 (SDC4); PDZ and LIM domain 2 (PDLIM2); talin 2 TLN2 (2 of 2)
Blood vessel development (GO:0001568: BP)
Forkhead box H1 (FOXH1); lysophosphatidic acid receptor-2(LPAR2 (2 of 2)); angiopoietin-2 (ANGPT2); Rho guanine nucleotide exchange factor-7(ARHGEF7(3 of 3)); melanoma cell adhesion molecule (MCAM1 of 2); angiopoietin-like-1(ANGPTL1 (2 of 2)); activin A receptor, typeI (ACVR1); activin A receptor type II-like-1(ACVRL1)
Proteinaceous extracellular matrix (GO:0005578: CC)
Netrin-4 (NTN4 (2 of 2)); sparc/osteonectin, cwcv and kazal-like domains proteoglycan (testican) (SPOCK2 (2 of 2)); collagens (COL9A2, COL5A2, COL11A1a, COL11A1b, COL27A1 (2 of 2)); ADAMTS8
Most genes identified as female-biased in the brain were expressed in both sexes but with significantly higher expression in females (Figure 3A). These transcripts were enriched with GO terms related to cell migration, cell adhesion, development, DNA replication, growth, glycolysis, and immune response (Table 3A, Additional file 9: Table S7). Several top female-biased transcripts encoded components of the proteinaceous extracellular matrix. For instance, genes encoding nidogens, laminins, fibronectins, collagens, as well as specific matrix metalloproteinases (Mmp-2-14) and members of disintegrin and metalloproteinase with thrombospondin motifs (Adamts) family were higher expressed in female brain (Figure 3B, Table 3A). Annotated genes with the strongest expression bias in the female brain included genes encoding peptide hormones- growth hormone-1 (Gh1), chorionic gonadotrophin beta 1 (Cgb1) and prolactin (Prl); and the calcium binding proteins parvalbumin-2 (Pvalb2) and calsequestrin-1 (Casq1a). Expression of the gene encoding teleost brain-specific aromatase, cytochrome P450 19A1b (Cyp19a1b), was 5-fold higher in the female than the male brain (Figure 3B, Additional file 6: Table S4).
We found similar numbers of male-biased and female-biased genes in the tail (Table 2, Figure 3C). GO terms related to signaling pathways, vesicle organization and transport, and transmembrane transport were enriched in the male-biased sequences (Table 3B, Additional file 10: Table S8). Several of the top male-biased genes encoded proteins with functions in pigment biosynthesis (Figure 3D, Additional file 7: Table S5, see below for more detail). Among female-biased genes, GO terms for cell-division, DNA replication, repair and recombination, glycolysis, and extracellular matrix components were enriched (Table 3B, Additional file 10: Table S8). Differentially expressed genes with growth-related functions included genes encoding mitotic cell-cycle factors - cyclin B1, cyclin A2, cyclin dependent kinase-1, and mini-chromosome maintenance (MCM) replication initiation factors (Figure 3D, Table 3B).
In gonads, 77% of all expressed genes showed sex-biased expression (Table 2, Figure 3E). We also found a number of genes with probable sex-limited expression in ovary or testis (Figure 3E). Male-limited and male-biased transcripts showed a greater magnitude of fold-change than the female-biased transcripts (Figure 3E, 3F). These included genes encoding some male-specific sex-development and differentiation associated proteins (e.g. DM-domain transcription factor Dmrt1, its paralog Dmrt2, and the 11-ketotestosterone biosynthesis enzyme Cyp11b2) (Additional file 8: Table S6); sperm associated antigens, ciliary and flagellar proteins (e.g., Spag17, Spag6, Tekt-1); spermatogenesis related - Spatc1l, Spata4; and testis expressed Tex9 (Figure 3F). Enriched GO-terms associated with male-biased genes included cilium assembly, spermatogenesis, and microtubule-based movement (Table 3C, Additional file 11: Table S9). Among the female-limited and female-biased sequences, we found genes encoding aromatase A (Cyp19a1a), the zona pellucida glycoproteins Zp1 and Zp2, oocyte specific proteins Zar1 and Zar1l (Figure 3F), ovarian folliculogenesis factors Gdf9 and Bmp15, and forkhead domain transcription factors Foxl2 and Foxr1 (Table 3C, Additional file 8: Table S9). Over-represented GO terms associated with female-biased genes were blood vessel development, regulation of BMP signaling pathway, amino acid transport, focal adhesion, cell migration involved in gastrulation, FGF receptor signaling, apical protein localization, regulation of body-fluid levels, and gas transport (Tables 3E, Additional file 11: Table S9).
A greater number of female-biased than male-biased genes showed a common direction of sex-bias in two or all three tissues (Table 2). Over-represented GO terms among genes with female-biased expression in both brain and tail included glycolysis, DNA replication and recombination as biological process terms, and proteinaceous extracellular matrix as cellular component term. While only a few genes showed male-biased expression in both brain and tail, we found that enriched terms related to transmembrane ion transport were common to both (Tables 3A, B, Additional file 12: Table S10).
Based on alignment positions of all genes on the currently available female draft genome sequence, we analyzed the distribution of all sex-biased genes (1.2 fold; FDR < 0.05) in comparison to all expressed genes (log2CPM > 2) along the guppy linkage groups. The total number of sex-biased genes per chromosome with their observed proportions and significance values for difference from expected proportions is described in Additional file 14: Table S12.
Linkage groups (LGs) with over-representation or under-representation of male-biased or female-biased genes
Male-biased genes over represented
Male-biased genes under represented
Female-biased genes over represented
Female-biased genes under represented
LG2, LG4, LG8, LG11, LG12, LG15, LG17, LG22
LG2, LG9, LG12, LG17, LG22
Comparison of d N /d S values for sex-biased genes to unbiased genes for brain, tail, and gonad tissues
0.0038 (< 0.0001)
0.00014 (< 0.0001)
Sex-linked genes may evolve at faster rates due to recombination differences between male and female germline, we therefore repeated the analysis using only autosomal genes and found similar rates of coding sequence evolution (Table 5). Also, magnitude and breadth of gene expression may be associated with functional constraints on coding sequence evolution ; we therefore further compared autosomal genes that were sex-biased in single tissues, multiple tissues, and had overall high expression (log2CPM > 5 i.e. 32 counts per million). Grouping by expression level or by expression breadth did not change the trend for higher d N /d S of genes with female-biased expression in brain or sex-biased expression in the gonads (Table 5).
Recent studies have shown that different assembly algorithms have varying strengths and limitations and a comprehensive assembly strategy should include the use of multiple assemblers [55–57]. While information from genomic coordinates assists in the assembly of full-length transcripts, at the same time genome-independent assemblies are free from biases caused by potential gaps and mal-orientations found in draft genomes. Therefore, we combined both assemblies to generate a non-redundant reference transcriptome composed of 74,000 CDS. Approximately, 24,000 CDS (~35%) were assigned as bona fide orthologs of published coding sequences. The remaining sequences may represent partially assembled sequences as well as incomplete CDS predictions and they may also include as yet unknown CDS (e.g., novel paralogs or splice variants), noncoding RNA, or assembly artifacts like chimeric transcripts. Our reference transcriptome provides a resource for investigating the genetics of complex adaptive traits such as guppy color patterns, life-history, and behavior [32, 38, 58].
Based on GO terms and orthology predictions, we attempted to relate our observations on sex-biased gene expression to sexually dimorphic phenotypic traits of the guppy. Pigment cells contributing to the adult male ornaments were expected to show a sex-biased expression mainly in adult skin, included as part of the tail in this analysis. Of the candidate genes associated with pigmentation, several were indeed higher expressed in male tails. A distinctive trait of the live-bearing female guppies is their lifelong growth, while male growth slows down after puberty [32, 59]. In concordance with this phenotypic difference, transcripts encoding cyclins and kinases characteristic of the mitotic phase of the cell cycle, DNA replication proteins and several growth factors were higher expressed in the female tail.
Female-biased expression of genes encoding cell-cycle and growth related hormones was also observed in the brain. Moreover, transcripts of the neurogenic zone associated aromatase, cyp19a1b, were highly expressed in the female brain but not the male brain, suggesting sexual dimorphism in adult neurogenesis in the guppy [60, 61]. We found a female bias in expression of many ECM components, which previously have been associated with neurogenesis and synaptic plasticity [62, 63]. Interestingly, greater brain plasticity in females as compared to male guppies has been postulated based on predator avoidance, kin-recognition, and mate choice differences in the wild [38, 42, 64, 65]. We also detected male-biased transcripts that encode neuropeptides and several transmembrane receptors in the brain, suggesting sex-differences in signal transduction. One example of such a male-biased transcript encodes the neuropeptide galanin, known to be involved in the neuroendocrine regulation of growth and reproduction in fish . Galanin neuropeptide and its receptor have also been shown to be more highly expressed in parts of the male versus female brain of sailfin mollies (Poecilia latipinna) [67, 68].
The highest degree of sex-bias in gene expression was found in the gonads, as expected in a gonochoristic organism. Expectedly, testis-biased transcripts encoded proteins with functions pertinent to testicular cells, e.g. spermatogenesis, sperm motility, and meiosis. Testis-specific and testis-biased expression of genes encoding Dmrt1 and Dmrt2, respectively, suggests a requirement of these transcription factors for testis maintenance . Ovary-biased transcription factors included the steroidogenesis regulators Foxl2 and Nr5a1 . Continuous FOXL2 activity is known to be required for suppression of trans-differentiation of ovarian cells into testicular cells in adult mice [71, 72]. The ovary-specific expression of the aromatase Cyp19a1a and testis-biased expression of Cyp11b2, which encodes for an enzyme for androgen 11-ketotestosterone biosynthesis , also indicates differences in sex-steroid synthesis in the two tissues. According to functional GO classification, female-biased genes were enriched for follicular vascularization factors, likely related to the lecithotrophic developmental strategy of the guppy . During oocyte growth in lecithotrophic species, the follicle is required for efficient transport of yolk precursors and probably amino acids and other metabolites from the blood to the maturing oocyte. After fertilization, the highly vascularized follicle persists as a placenta that is essential for osmoregulation, gas exchange, waste disposal, and transport of some essential factors [74, 75].
We found a significant enrichment of female-biased genes on the putative X-chromosome, LG12. Sex-biased genes have been reported to accumulate on differentiated sex-chromosomes of many species. Enrichment of female-biased genes on X-chromosomes has been reported in species with heterogamous males, e.g. several Drosophila species [21, 76], mouse , and the nascent X-chromosome of the stickleback (Gasterosteus aculeatus) . Similarly, enrichment of male-biased genes has been found on the Z-chromosomes of species with heterogamous females, e.g. zebra finch (Taeniopygia guttata) and chicken . In guppies, the majority of the sex chromosome is pseudo-autosomal, yet the X and Y chromosome show genetic and cytological distinctions [54, 78]. Although differentiation between X- and Y-chromosomes is not comparable to the situation in mammals, birds, or drosophilids, the over-representation of ovary-biased genes and under-representation of testis-biased genes on the guppy LG12 indicates sex-specific selection pressures even in the absence of a truly hemizygous state in males. Previous studies have indicated reduced synaptic pairing  and reduced recombination  between X- and Y-chromosomes during male gametogenesis in guppies. This may lead to accumulation of deleterious mutations in Y-linked alleles of genes on LG12. Ovary-biased or ovary-specific genes are likely not needed in males and therefore mutations in these genes will persist on the Y-chromosome, while mutations in testis-biased genes and other non-biased genes will be selected against. Further analysis of recombination frequencies and gene order along the length of the sex chromosomes, coupled with comparisons across multiple populations, will enable better understanding of the effect of genomic location of sex-biased genes.
Our comparisons of d N /d S of sex-biased and unbiased genes in the guppy revealed elevated coding sequence change of testis and ovary biased genes, and female-biased genes expressed in the brain and those co-expressed with a female-bias in the brain and at least one of the other tissues. Current knowledge on protein evolution suggests that sex-biased genes in reproductive and non-reproductive tissues show accelerated evolution in many species including Drosophila, mouse, and chicken [18, 19, 79, 80]. Sex-biased genes may show rapid divergence due to their evolution under sexual selection. Additionally, accelerated sequence divergence may also occur under sex-specific natural selection, or relaxed purifying selection on genes that have reduced functional pleiotropy [15, 81]. Our results indicating rapid evolution of sex-biased genes in gonads driven by testis-biased genes are in concordance with the results obtained from other vertebrates [18, 82, 83]. Rapid divergence of reproductive proteins driven by testis-specific proteins may be involved in sexually-antagonistic selection , post-copulatory sexual selection resulting from co-adaptation , or kinship recognition and incipient speciation in guppies . These processes may be important in guppies given their highly promiscuous mating system with coercive internal fertilization by males and long-term sperm storage in females .
We also found a higher d N /d S ratio in female-biased genes expressed in the brain. While considerable evidence suggests that sexual selection in guppies is driven by the multivariate mating preferences of females, male-male competition, male mating tactics, and male mate choice may also be under selection [37, 38, 87, 88]. An association between molecular evolution of female-biased genes and sexual selection on behavioral traits was not easy to elucidate as these genes were expressed in both sexes and co-expressed in multiple tissues, suggesting some pleiotropy in function. Rapid evolution of female-biased genes with growth-related and metabolic functions may be pertinent to the sexual size dimorphism seen in this species and may be driven by natural selection on life-history traits .
Unexpectedly, we also found signatures of rapid evolution in female-biased genes whose expression was not tissue-specific. Usually a broad expression profile indicates pleiotropic functions that would constrain divergence of coding genes [19, 90]. This prediction is, however, not necessarily cogent for fish, where the teleost-specific whole genome duplication allows for evolving sub-functionalization or even some redundancy when co-orthologs are preserved . Furthermore, many of the co-expressed female-biased genes identified in our study encode ECM components, cell-cycle factors, and glycolytic enzymes. These proteins have conserved functions across all tissues and therefore may not be pleiotropic. Moreover most of these proteins are located in the cytosol or in the extracellular region where adaptive evolution has been described to be more relaxed [90, 91]. Conversely, no difference in evolutionary rate was found between male-biased genes co-expressed in the brain and tail and unbiased genes. Many of these genes encoded neuropeptides, transmembrane receptors, and gated ion-channels that evolve under structural and functional constraints of ligand-receptor specificity and transmembrane localization . Therefore, these proteins are likely to have low tolerance of mutations and more conserved evolution.
Our analyses of sex-biased gene expression in guppies revealed differences that are likely to be pertinent for the mechanisms underlying its sexual dimorphism. The observed female-biased expression of growth-related genes in the three tissues investigated could reflect the life-long growth observed in female guppies. At the same time, elevated male-biased expression of genes known to be relevant for pigment cell differentiation was limited to the tail, the tissue including part of the adult skin. As expected, sperm-specific and cell-cycle-relevant transcripts were highly expressed in the testis, and the expression profile of the ovary signifies maternal provisioning of this lecithotrophic species. Correlations between gene expression and phenotypic traits will remain speculative in guppies until methods of experimental gene gain and loss of function can eventually be established in this system.
We detected accumulation of ovary-biased genes on the putative X-chromosome, supporting the hypothesis that genes advantageous to only one sex accumulate on the differentiating sex chromosomes. We also observed more rapid evolution of testis-biased genes, possibly indicating increased sexual selection on males. However, the observed rapid evolution of genes with female-biased expression in brain and other tissues not seen in males may be confounded by differences in the biological functions and cellular locations of male- and female-biased genes. It is probable then, that there are differences in selection on protein sequences of males and females independent of the breadth of tissue expression.
This study was carried out in accordance with the German Protection of Animals Act (§ 11 Abs. 1 Nr. 1 a und b TierSchG). All animal experiments were permitted by the Regierungspräsidium Tübingen (approval ID 35/9185.46). Sample tissues were prepared from laboratory-reared guppies that were descendants of wild fish caught in 2003 from a low-predation population in Quare river, East Trinidad [45, 92]. The fish were reared under uniform conditions of food, water, light, and density. Mature adult guppies between 5–6 months of age were isolated and kept without food in fungicide treated water for 44–48 hours prior to dissections. Fish were euthanized using 0.1% (w/v) tricaine (ethyl 3-aminobenzoate methanesulfonate salt) solution pH 7. Brain, eyes, liver, spleen, skin, tail (the post-anal tissue up to the beginning of the tail fin, containing adult skin, skeletal muscle, dorsal cord, bone, and cartilage), and gonads were isolated from euthanized adult males and females. Female tissues were prepared from virgin adults that were separated from males at the age of 3–4 weeks to avoid pregnancy and sperm storage. Whole embryos at late-eyed to very late-eyed stages of development  were isolated from gravid females. A small fin-clip was taken from each embryo for genotyping the sexes. All samples were washed with cold PBS (kept at 4°C), then frozen in liquid nitrogen, and stored at −80°C.
Four Illumina cDNA libraries were prepared separately from female and male late-eyed stage embryos, and adult female and male tissues (brain, eyes, liver, spleen, skin, tail, and gonad tissues). Embryos were first genotyped using genomic DNA isolated from fin-clips and markers 229 and 230 with sex-specific single nucleotide polymorphisms (SNPs) in Quare population .
All tissue samples were homogenized in TRIzol® reagent (Invitrogen) using a Polytron® homogenizer (PT 1200, Kinematica AG, Switzerland) and total RNA was extracted from the Trizol homogenate according to manufacturer’s instructions. After removal of contaminant DNA, using DNaseI (Invitrogen), purified RNA was quality-checked and quantified (Nanodrop ND-2000, ThermoScientific peqlab®). For male and female adult libraries, a total of 75 μg RNA was prepared by pooling 15 μg of RNA isolated from each tissue. For male and female embryo libraries, 75 μg total RNA was isolated from 15 individual embryos each. Then purified polyA + mRNA (Dynabeads® Oligo(dT), Invitrogen) was used for preparation of paired-end RNA libraries with 200-300 bp insert size, using the mRNA-seq Sample Preparation Kit (Illumina, San Diego, CA) or the NEBNext® mRNA Library Prep Reagent Set for Illumina (NEB), according to manufacturer’s instructions. Library quality and concentrations were assessed using the Agilent DNA 1000 Bioanalyser Assay (Agilent Technologies, Germany). Each library was sequenced on a separate GAIIx lane (Illumina, San Diego, CA, read length 101 bp). These four datasets are referred to as female and male adult (Fadult, Madult) and female and male embryo (Fembryo, Membryo; Figure 1).
Barcoded cDNA libraries were prepared for quantitative analysis of gene expression differences. Tissues from adult male and female brain and eyes, tail, and gonads (ovaries from virgin females or testes from males) were isolated as indicated in Figure 2. All tissues were individually homogenized using steel beads for disruption in TRIzol Reagent (Invitrogen, Carlsbad, CA, USA). Total RNA was extracted from the TRIzol homogenate using DirectZol RNA extraction kits with in-column DNaseI treatment. Purified total RNA was quality-checked on agarose gels and quantified using the Qubit RNA Assay Kit (Invitrogen, Carlsbad, CA, USA). Six biological replicates were prepared for each tissue and sex, except the female brain, which consisted of 7 biological replicates and two technical replicates (but only 6 biological replicates were used for quantitative comparisons). All samples were randomized and individually barcoded during library preparation using TruSeq mRNA-seq Sample Prep Kit (Illumina, San Diego, CA, mRNA-seq Sample Prep Manual v2 protocol). In total, 39 paired-end libraries were prepared and sequenced on 3 lanes of the HiSeq™ 2000 (Illumina, San Diego, CA, read length 101 bp, 13 libraries per lane). The barcoded cDNA libraries from adult tissues are referred to as: Fbrain, Mbrain, Ftail, Mtail, Fgonad, Mgonad (Figures 1 and 2).
The types of tissues, number of individuals, types of libraries and sequence datasets are summarized in Additional file 1: Table S1.
The resulting reads in the non-barcoded datasets were filtered for low complexity using SHORE v0.6  and for PCR duplicates by removing read pairs that matched 60 bp of both reads of a pair and keeping unique pairs and 3 duplicates with highest quality scores (customized perl script). We trimmed homopolymer sequences (polyA/T/G/C) over 22 bp length using Cutadapt v1.2.1  and low-quality nucleotides using CON DE TRI v2.2  with a phred20 quality cutoff, 35 bp length cutoff and other default parameters. The barcoded dataset used in expression analysis was filtered only for low quality but trimmed similarly.
A genome-independent transcriptome was assembled using Trinity (trinityrnaseq_r2012-06-08)  with minimum k-mer coverage of 2 and default parameters. In order to maximize the k-mer overlap and to achieve high coverage for rare transcripts, we combined all datasets (Fadult + Madult + Fembryo + Membryo + Fbrain + Mbrain + Ftail + Mtail + Fgonad + Mgonad) and assembled a single de novo reference trancriptome with Trinity (Figure 1).
A genome-guided transcriptome assembly was compiled using TOPHAT – CUFFLINKS – CUFFMERGE v2.0.4 [98, 99] with default parameters using a draft assembly of the female guppy genome (Künstner et al., in preparation). Reads from each RNA-seq sample were first individually mapped to the reference genome using TOPHAT2 to retain tissue-specific splicing information. The resulting alignment files were analyzed by Cufflinks to generate a transcriptome assembly for each dataset (Fadult, Madult, Fembryo, Membryo, Fbrain, Mbrain, Ftail, Mtail, Fgonad, and Mgonad). These assemblies were then merged to give a combined assembly with CUFFMERGE (Figure 1).
Open reading frames (ORFs) were predicted using the program Trans Decoder from Trinity. Predicted coding sequences (CDS) over 150 bp long were clustered if they were greater than 90% similar and the longest sequence was kept in a non-redundant database using Cd-Hit-Est v4.6 [100, 101].
Orthologous genes to other vertebrate species were identified using translated CDS for both genome-guided and genome-independent assemblies. Sequences from Danio rerio (zebrafish), Gasterosteus aculeatus (stickleback), Oryzias latipes (medaka), Xiphophorus maculatus (platyfish), Oreochromis niloticus (tilapia), Takifugu rubripes (fugu), Tetraodon nigriviridus (tetraodon), Gadus morhua (cod), Homo sapiens (human), and Mus musculus (mouse) were downloaded from Ensembl (Release 71) . Single-copy (1:1) orthologs were identified using Protein Ortho v4.26  (parameters: NCBI Blastp] v2.2.21, E-value <1 × 10−10, alignment connectivity: 0.8, coverage: 40%, identity: 30%, adaptive similarity: 0.95, including pairs: 1).
Pooled paired-end reads from all sequenced datasets were normalized using Diginorm with default parameters for single-pass normalization. Reads from the normalized dataset were aligned to the genome-guided and the genome-independent assemblies using BOWTIE2 v2.0.04  (default parameters for sensitive local alignment).
We merged the genome-independent and genome-guided assemblies by pooling the predicted CDS from both assemblies followed by clustering sequences with greater than 90% identity using Cd-Hit-Est to create a guppy reference transcriptome (GRT).
Annotations were found using NCBI BlastX v2.2.25 and the NCBI non-redundant protein database . Functional categories were assigned by mapping GO terms using Blast2GO® v2.7.0 . For simplicity we refer to the predicted CDS of the guppy reference transcriptome as genes in the text.
Genomic coordinates of genes in the reference transcriptome were obtained by aligning them against the repeat-masked draft female genome using GMAP v2012-07-20 . In the case of ambiguous alignments, the alignments with the highest total coverage and identity were kept (total 607).
Each barcoded sequenced library from the organ datasets (Fbrain, Mbrain, Ftail, Mtail, Fgonad, Mgonad) was individually aligned to genes of the guppy reference transcriptome using BOWTIE2 v2.0.04. Mapped reads were counted using eXpress v1.3.1 . Read counts from six individually barcoded biological replicates per tissue were used for differential expression analysis between male and female tissues using the Bionconductor package edgeR v3.0.8 . First low abundance genes with less than two counts per million mapped reads (<2 CPM/sample) across six samples were removed. Read counts were normalized for sequencing depth using TMM normalization  and differential expression between sexes was tested with a modified exact test implemented in edgeR  and corrected for multiple testing. Genes with significant expression difference between the sexes (FDR < 0.1 or if mentioned FDR < 0.05) were classified as sex-biased and those without a significant difference as ‘unbiased’. Using an FDR cut-off of 0.1 sex-biased sequences showed at least a 1.2 fold difference (log2FC > 0.3 or < −0.3) in expression between the sexes. Genes with sex-specific functions may have varying levels of expression divergence in different tissues [12, 19]. Therefore, we further categorized the sex-biased genes by fold-change, keeping genes with greater than median-fold difference in expression between sexes (log2FC: Male/Female) in each study tissue. These median-fold cutoffs were: 1.5 fold in the brain (log2FC > 0.6 or < −0.6); 1.7 fold in the tail (log2FC > 0.8 or < −0.8); and 3.2 fold in the gonad (log2FC > 1.8 or < −1.8). Enrichment of GO terms between sex-biased and unbiased genes per tissue was determined using a Fisher’s exact test with the Elim algorithm (p < 0.01, and number of sequences > 3) in the R package: topGO v2.10.0 .
Non-random chromosomal distribution of male- or female-biased genes expressed in a tissue was tested with a χ2-test. P-values were corrected for multiple testing using the Benjamini-Hochberg method . The expected distribution was calculated by assuming that sex-biased genes are randomly distributed across chromosomes and that their representation on a particular chromosome is proportional to the number of expressed genes on that chromosome. In the brain the average number of male-biased genes was significantly lower than the average number of female-biased ones; therefore for this tissue we calculated the expected frequency of male- and female-biased genes using their respective averages. We conducted this analysis twice, (1) where sex-bias was defined as a greater than the 1.2-fold-change between the sexes (FDR < 0.05) and (2) where sex-bias was defined as greater than median fold difference (FDR < 0.1). All comparisons were tested using statistical tests implemented in R package Stats version 2.15.2 .
Orthologous amino acid sequences between stickleback, medaka, and guppy obtained from Protein Ortho were aligned using Mafft v7.017b  and back-translated to nucleotide sequences for subsequent analysis. All alignments are available upon request. Substitution rates were estimated separately for synonymous (d S ) and non-synonymous (d N ) substitutions using a maximum likelihood method implemented in the Codeml program (model = 1, user tree specified according to the phylogeny) in the Paml package v4.1 . We excluded all alignments shorter than 150 bp or with d S larger than 2 to minimize statistical artifacts from short sequences and saturation effects in d S . Mean values of d N /d S of male-biased and female-biased genes were compared to unbiased genes with significant expression (log2CPM > 2) per tissue.
All comparisons were tested using the non-parametric Mann–Whitney test as well as permutation tests with 1000 replicates (data not shown) using R version 2.15.2 .
The RNA-seq reads reported in this study have been deposited in the National Center for Biotechnology Information Short Reads Archive, http://www.ncbi.nlm.nih.gov/sra (Study Accession ID: SRP033586). The predicted CDS of the guppy reference transcriptome are available from our institute’s website: ftp://ftp.tuebingen.mpg.de/ebio/publication_data/esharma/guppy_trans/trin_cuff_v14_cdhit90.fa.gz
Expressed sequence tags
Quantitative trait loci
Single nucleotide polymorphism
Open reading frame
Guppy reference transcriptome
Counts per million
False discovery rate
Fragments per kilo base per million
Number of synonymous substitutions per synonymous site
Number of non-synonymous substitutions per non-synonymous site
Sex determining locus
NCBI non-redundant protein database.
We thank David Reznick for guppy specimens. We thank Margarete Hoffmann, Christa Lanz, and Jens Riexinger for help in library preparation and sequencing. We thank Stefanija Topuz for fish care, Dino Jolic for help in analysis and Gunnar Rätsch and Daniel Koenig for suggestions and discussions. We thank Rebecca Schwab for comments and help with writing.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.