De Novo Assembly of the Manila Clam Ruditapes philippinarum Transcriptome Provides New Insights into Expression Bias, Mitochondrial Doubly Uniparental Inheritance and Sex Determination Fabrizio Ghiselli, 1, * Liliana Milani, 1 Peter L. Chang, 2 Dennis Hedgecock, 3 Jonathan P. Davis, 4 Sergey V. Nuzhdin, 2 and Marco Passamonti 1 1 Dipartimento di Biologia Evoluzionistica Sperimentale, Universita ` di Bologna, Bologna, Italy 2 Program in Molecular and Computational Biology, Department of Biological Sciences, University of Southern California 3 Marine and Environmental Biology Section, Department of Biological Sciences, University of Southern California 4 Taylor Shellfish Farms, Quilcene, Washington * Corresponding author: E-mail: [email protected] Associate editor: Richard Thomas Abstract Males and females share the same genome, thus, phenotypic divergence requires differential gene expression and sex- specific regulation. Accordingly, the analysis of expression patterns is pivotal to the understanding of sex determination mechanisms. Many bivalves are stable gonochoric species, but the mechanism of gonad sexualization and the genes involved are still unknown. Moreover, during the period of sexual rest, a gonad is not present and sex cannot be determined. A mechanism associated with germ line differentiation in some bivalves, including the Manila clam Ruditapes philippinarum , is the doubly uniparental inheritance (DUI) of mitochondria, a variation of strict maternal inheritance. Two mitochondrial lineages are present, one transmitted through eggs and the other through sperm, as well as a mother- dependent sex bias of the progeny. We produced a de novo annotation of 17,186 transcripts from R. philippinarum and compared the transcriptomes of males and females and identified 1,575 genes with strong sex-specific expression and 166 sex-specific single nucleotide polymorphisms, obtaining preliminary information about genes that could be involved in sex determination. Then we compared the transcriptomes between a family producing predominantly females and a family producing predominantly males to identify candidate genes involved in regulation of sex-specific aspects of DUI system, finding a relationship between sex bias and differential expression of several ubiquitination genes. In mammalian embryos, sperm mitochondria are degraded by ubiquitination. A modification of this mechanism is hypothesized to be responsible for the retention of sperm mitochondria in male embryos of DUI species. Ubiquitination can additionally regulate gene expression, playing a role in sex determination of several animals. These data enable us to develop a model that incorporates both the DUI literature and our new findings. Key words: Ruditapes philippinarum , de novo, transcriptome, doubly uniparental inheritance, sex bias, sex determination. Introduction Males and females undergo different selective pressures, some operating in opposite directions. Because both sexes share the same genome (except for sex chromosomes, where present), phenotypic divergence requires sex- specific regulation ( Ellegren and Parsch 2007 ; Arnold et al. 2009 ), and males and females, even with the same set of genes, show differences in gene expression or use al- ternative splice forms ( Long et al. 1995 ; Nuzhdin et al. 1997 ; Jin et al. 2001 ; McIntyre et al. 2006 ; Foley et al. 2007 ; Chang et al. 2011 ). Overall, sex-related differences in gene expres- sion were observed across a wide range of taxa ( Ellegren and Parsch 2007 ). For example, over 12% of the germ line transcripts of Caenorhabditis elegans showed a sex bias, and expression analyses on whole Drosophila melanogaster body showed that the proportion of genes presenting a sex bias is around 57% ( Jin et al. 2001 ; Arbeitman et al. Sex Determination in Bivalves and DUI Many bivalves are stable gonochoric species, but the mech- anism of gonad sexualization and the genes involved are still unknown ( Paz et al. 2005 ; Breton et al. 2007 ). During the period of sexual rest, a gonad is not present and sex cannot be determined. Every year in the reproductive sea- son, testis and ovary develop from a group of germ cells ( Devauchelle 1990 ; Milani L, Ghiselli F, Maurizii MG, and Passamonti M, in preparation), and sex can be determined by detecting sperm or oocytes microscopically. In addition, heteromorphic sex chromosomes appear to be absent from bivalves ( Sastry 1979 ; Borsa and Thiriot-Quie ́ vreux 1990 ). DUI is a mechanism associate d with germ line differenti- ation in some bivalves, including R. philippinarum ,anotewor- thy variation of the strict maternal inheritance. In DUI species, two mitochondrial lineages ar e present, one transmitted through eggs (called F, for female-transmitted) and the other transmitted through sperm (call ed M, for male-transmitted). In DUI, both sexes inherit F mitochondria from the mother, whereas M mitochondria are transmitted from father to sons only ( Breton et al. 2007 ; Passamonti and Ghiselli 2009 ). Thus, two different mitochondrial genomes, with an unexpectedly high level of sequence divergence, up to 52% ( Doucet- Beaupre ́ et al. 2010 ), are detectable. In addition, a high var- iability of progeny sex ratio was observed in Mytilus DUI species: Some females produce female-biased offspring, others male-biased and still others a 1:1 ratio ( Saavedra et al. 1997 ; Kenchington et al. 2002 , 2009 ; Cogswell et al. 2006 ). The ex- istence of lineages presenting skewed sex ratios in DUI ani- mals has been proposed to be a peculiarity of their sex determination mechanism. Exp erimental evidence suggests that control over sex ratio is exercised by the mother’s nuclear genome ( Saavedra et al. 1997 ; Kenchington et al. 2002 , 2009 ; Cogswell et al. 2006 ). Indeed, matings of the same female with different males always give the same sex ratio, but matings of the same male with different females result in different sex ratios. The genetic factors involved are so far unknown, but probably sex determination in bivalves is oligogenic, with multiple coexisting genes ( Kenchington et al. 2002 ). In embryos of analyzed DUI species, sperm mitochondria follow two different distribution patterns. In males, they aggregate near the cleavage furrow at the first cell division, eventually segregating into the primordial germ cells, whereas in females they are dispersed among blastomeres and degraded ( Cao et al. 2004 ; Obata and Komaru 2005 ; Cogswell et al. 2006 ; Kenchington et al. 2009 ; Milani L, Ghiselli F, and Passamonti M, in preparation). Sperm mito- chondria degradation in mammal embryos is mediated by ubiquitination ( Sutovsky et al. 2000 ; Here, we report the first whole transcriptome analysis by RNA-Seq performed to identify genes involved in bivalve sex determination and DUI. The characterization of genes associated with reproduction and the analysis of their ex- pression pattern and polymorphism can provide insight in- to molecular mechanisms regulating sex determination. We compared the transcriptomes of males and females and identified 1,575 genes with strong sex-specific expres- sion and 166 sex-specific single nucleotide polymorphisms (SNPs), obtaining preliminary information about genes that could be involved in sex determination. Furthermore, for the first time in a DUI species outside the Mytilus complex, we confirmed the presence of sex-biased families. Then, we compared the transcriptomes between a family producing predominantly females and a family producing predomi- nantly males to identify candidate genes involved in the regulation of sex-specific aspects of DUI system. Finally, we produced a model that is consistent with the DUI lit- erature and with our new transcriptomic data. Materials and Methods Clam Families In 2006, Taylor Shellfish Farm (Quilcene, WA) generated 29 families of clams by pairwise crosses of animals originally caught in the wild. In the summer of 2009, the sex ratio of these families was determined by relaxing the clams with MgSO 4 and taking a needle biopsy of gonad tissue for mi- croscopic examination ( table 1 ). Additional individuals from two families, which were selected for their strong sex bias, 023 (female-biased) and 025 (male-biased), were shipped to the University of Southern California, where the live animals were opened and sexed. The bodies of three males and three females from each family were frozen in liquid nitrogen for eventual preparation of cDNA libraries. Library Preparation RNA purification, cDNA synthesis, and Illumina library con- struction were performed using the protocols of Mortazavi et al. (2008) , with the following modifications. Total RNA, mRNA, and DNA were quantified using a Qubit fluoro- meter (Invitrogen). mRNA fragmentation was performed using Fragmentation Reagent (Ambion) for a 3 min and 50 s incubation at 70 ° C and subsequently cleaned through an RNA cleanup kit (Zymo Research). Additional DNA and gel purification steps were conducted using Clean and Con- centrator kits (Zymo Research). Each sample individual was barcoded following the Illumina protocol. Two technical replicates were generated per individual for paired-end 71-bp reads on an Illumina Genome Analyzer II, producing a total of 2 technical replicates
3 biological replicates
2 sexes
2 families 5 24 samples. Table 1. Family Sex Ratios. Family ID Total Numbers Males Females Number of Sexed Percentage of Sexed Percentage of Males 032 19 1 11 12 63 8 023 a 60 8 38 46 77 17 014 18 4 9 13 72 31 001 17 5 10 15 Short-Read Sequencing and De Novo Assembly Across all samples, 90 million (M) paired-end reads, 71 ba- ses long, were obtained. All raw Illumina FastQ sequences are available for download at the NCBI Short Read Archive under the accession number SRA037984.1. Because de novo assembly of transcriptomes from nonmodel species lacking genome sequences can be sensitive to sequencing errors, it is critical that the reads used to generate contigs have the highest sequencing quality. Reads were removed from consideration in the de novo assembly if the read had a terminal ‘‘phred’’ ( Ewing and Green 1998 ) quality value less than 15, or if the read contained more than two un- known nucleotides (i.e., N ). Reads were also filtered if they were similar to known polymerase chain reaction primer and Illumina adapter sequences. Since the M and F mito- chondrial genomes are available (GenBank accession num- bers AB065374 and AB065375), additional reads were removed when aligned with six or fewer mismatches to these sequences. After these four filtering steps, 32.9 M paired-end reads and 8.1 M single-end reads were kept for further analyses ( supplementary table S1 , Supplemen- tary Material online). The 41 M retained reads were assembled with Velvet (version 1.0.15) ( Zerbino and Birney 2008 ), in conjunction with a custom post-processing algorithm capable of retain- ing information from alternative splices. Velvet was run un- der the following settings with a kmer length of 35: -cov_cutoff auto -max_branch_length 0 -max_divergence 0 -max_gap_count 0 -read_trkg yes. Sequenced reads that were kept as pairs and not filtered out together or sepa- rately were treated as ‘‘-shortPaired’’ with insert length of 105 bases and standard deviation of 40 bases. Single- end reads that were not filtered out were treated as ‘‘-short.’’ ASplice was run with default parameters. These algorithms utilize de Bruijn graphs to assemble short reads into contigs, using sequence overlap information, until the contigs can no longer be extended. The assembly resulted into 35,784 contigs clustering within 22,886 nodes. These nodes represent genes and their isoforms identified in the assembly. The assembled sequences are available on the NCBI Transcription Shotgun Assembly Database (Gen- Bank Accession Numbers: JO101212-JO124029). With the set of de novo sequences serving as a reference, reads from each of the individual samples were mapped using the Burrows-Wheeler Alignment tool ( Li and Durbin 2009 ). The number of reads that mapped to each gene was tabulated and normalized to calculate fragments per kilo- base of exon per million fragment sequenced (FPKM). Ad- ditional normalization among all samples was performed, using the Trimmed Mean of M values protocol outlined in Robinson and Oshlack (2010) , which takes into account differences in overall RNA populations across biological samples and is one of several methods used to evaluate RNA-sequencing data. Normalization was implemented us- ing the edgeR package in R ( Robinson and Smyth 2007 ). To detect differential expression, we used FPKM values as the dependent variable in gene-specific mixed linear models implemented in R. Sex, family, and sex-by-family interac- tion were fixed effects, whereas replicates were random ef- fects. The false discovery rate was used to account for multiple testing ( Benjamini and Hochberg 1995 ); a cutoff of 0.05 was applied to call effects significant. The signifi- cance of differential gene enrichment between groups was tested using a Wilcoxon rank sum test. All statistical analyses and graphs evaluating consistency between replicates and genotypes were produced using R v2.13.0. Gene Functional Annotation and Classification Blast2GO v.2 ( Go ̈tz et al. 2008 ) and WEGO ( Ye et al. 2006 ) were used to obtain dene ontology (GO) annotations. Genes were also annotated using a BLASTX ( Altschul et al. 1990 ) search to the nonredundant GenBank CDS transla- tions þ PDB þ SwissProt mussel are available, but the divergence time with R. phil- ippinarum was estimated between 542 and 488 Ma ( Plazzi and Passamonti 2010 ). For this reason, we allowed a higher flexibility and chose the annotation with the highest BLAST score as long as the span of the alignment was greater than 80% of the length of the gene under query. For genes that did not report any hits, we lowered the minimum span to 40% of the length, choosing the annotation with the high- est BLAST score, having Expected value , 1.00
10  5 . The GOstat package ( Beissbarth and Speed 2004 ) was used to identify overrepresented GO categories in groups of tran- scripts ( P , 0.01). InterProScan version 4.8 ( Hunter et al. 2009 ) was used to identify functional conserved domains of reproductive and ubiquitination genes. Sequence Polymorphism Analysis Representative transcript sequences were identified using a global multiple sequence alignment of all contig sequen- ces for each node. For each sample, SNPs were identified with reference to the de novo assembled reference sequence, using SAMtools ( Li et al. 2009 ). Given the nature of the assembly, the SNP data were calculated in a conservative and parsimonious way: Sites with less than 5
coverage were discarded, positions with a phred score lower than 15 were excluded, and indels were not taken into account. All assembled sequences were then aligned and analyzed with the VariScan 2.0 software ( Hutter et al. 2006 ) in order to compute polymorphism data. A block data file was generated to specify gene boundaries in the alignment in order to calculate the statistics for each gene. The program was run under the following settings: RunMode 5 12, UseMuts 5 1, CompleteDeletion 5 0, Fix- Num 5 1, NumNuc 5 9, SlidingWindow 5 0. This config- uration reported the number of segregating sites ( S ), total number of mutations ( g ), the number of singletons, nucle- otide diversity ( p ), Watterson’s estimator of nucleotide di- versity per site ( h ), Tajima’s D statistic, Fu & Li’s D* and F* . Sequence polymorphism was analyzed between the follow- ing differentially expressed gene categories: family-biased genes versus family-unbiased genes, sex-biased genes versus sex-unbiased genes, male-biased genes versus female-biased genes, and reproductive genes versus male-biased genes. The reproductive gene group included genes annotated by Blast2GO under the ‘‘Reproduction’’ category. We used R Table 3. Mean Expression of Genes. Sex-Biased Genes ( 6 SE) Male-Biased