Bivalves exhibit an astonishing diversity of sexual systems and sex-determining mechanisms. They can be gonochoric, hermaphroditic or androgenetic, with both genetic and environmental factors known to determine or inﬂuence sex. One unique sex-determining system involving the mitochondrial genome has also been hypothesized to exist in bivalves with doubly uniparental inheritance (DUI) of mtDNA. However, the link between DUI and sex determination remains obscure. In this study, we performed a comparative gonad transcriptomics analysis for two DUI-possessing freshwater mussel species to better understand the mechanisms underlying sex determination and DUI in these bivalves. We used a BLAST reciprocal analysis to identify orthologs between Venustaconcha ellipsiformis and Utterbackia peninsularis and compared our results with previously published sex-speciﬁc bivalve transcriptomes to identify conserved sex-determining genes. We also compared our data with other DUI species to identify candidate genes possibly involved in the regulation of DUI. A total of12,000 orthologous relationships were found, with 2,583 genes differentially expressed in both species. Among these genes, key sex-determining factors previously reported in vertebrates and in bivalves (e.g., Sry, Dmrt1, Foxl2) were identiﬁed, suggesting that some steps of the sex-determination pathway may be deeply conserved in meta- zoans. Our results also support the hypothesis that a modiﬁed ubiquitination mechanism could be responsible for the retention of the paternal mtDNA in male bivalves, and revealed that DNA methylation could also be involved in the regulation of DUI. Globally, our results suggest that sets of genes associated with sex determination and DUI are similar in distantly-related DUI species. Key words: sex determination, mitochondrial DNA, doubly uniparental inheritance, comparative transcriptomics, Bivalvia. Introduction determine sex in most bivalve species studied to date Bivalves show extreme diversity in their reproductive systems, (reviewed in Coe 1943; Chavez-Villalba et al. 2011; Breton ranging from functional (simultaneous) hermaphroditism, al- et al. 2017). However, because sex determination has been ternative sexuality (sequential hermaphroditism), to strict gon- studied in great detail only in oysters, the genetic and envi- ochorism (i.e., species that exist as separate males and ronmental factors that control sexual diversiﬁcation in bivalves females). Both genetic and environmental factors appear to are still poorly known. An important reason for studying sex The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact firstname.lastname@example.org Genome Biol. Evol. 10(2):577–590. doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 577 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE determination in bivalves relates to their biology and ecology: orf in the F mtDNA and M-orf in the M mtDNA, both of which with 25,000 living species (www.bivatol.org), they are suf- are relatively conserved across species within a family (Breton ﬁciently diverse to provide a rich source of material to better et al. 2009, 2011; Milani et al. 2013; Mitchell et al. 2016). understand the evolution of sex and sex determination in These genes were hypothesized to be key elements of a sex- general, and to provide unique examples of sex-determining determination system involving mitochondria, that is, they mechanisms, including the only possible example of a sex- may play a role in the maintenance of a gonochoric repro- determining system involving the mitochondrial genome in ductive system with separate male and female sexes (Breton animals (e.g., Breton et al. 2011). et al. 2011). This hypothesis came from the observation that A variety of approaches have been undertaken to elucidate gonochorism in freshwater unionoid mussels is absolutely cor- the mode of sex determination in bivalves. Globally, these related with the presence of DUI and of these novel sex- studies demonstrated an absence of heteromorphic sex chro- speciﬁc proteins, whereas closely related hermaphroditic spe- mosomes, an effect of several environmental factors on sex- cies lack the M mtDNA (i.e., possess strict maternal inheri- ratio (i.e., temperature, food availability, and pollutants, such tance of mtDNA) and have macromutations in the F-orf as heavy metals organochlorines and exogenous steroids), gene in their F mtDNA (Breton et al. 2011). If true, this would and a polygenic architecture of sex determination (Breton make DUI the ﬁrst animal sex-determination system involving et al. 2017). Genes homologous to sex-determining pathway the mitochondrial genome, and would explain its long-term genes in other animals have been identiﬁed in several bivalve persistence in bivalves. However, this hypothesis still has not species. For example, homologs of Sry (sex determining been rigorously tested. region-y)-box 30 (Sox 30)and Dmrt1 (doublesex and mab-3 Ghiselli et al. (2012) published the ﬁrst whole transcrip- related transcription factor 1), known for their role in male sex tome analysis by RNA-Seq to better understand the mecha- determination in nematodes, fruit ﬂies and vertebrates (Wallis nisms underlying DUI and sex determination in the Manila et al. 2008; Kopp 2012), have been found in the marine clam clam R. philippinarum. Because previous studies demon- Ruditapes philippinarum, and in some species of scallops and strated that elimination of paternal mitochondria in mammal oysters (Yu et al. 2011; Ghiselli et al. 2012; Llera-Herrera et al. species with strict maternal inheritance of mtDNA is depen- 2013; Teaniniuraitemoana et al. 2014; Zhang et al. 2014; Li dent on a mechanism involving the ubiquitin–proteasome sys- et al. 2016). Several other genes known to act in sex deter- tem, the principal system for protein degradation (Sutovsky mination in other animals or during early gonadal differenti- et al. 1999), it was proposed that a modiﬁcation of the ubiq- ation have been identiﬁed in bivalves, such as Oyvlg and uitination process in male bivalves with DUI would allow vasph, that is, homologues of vasa, which is a gene involved sperm mitochondria and their M genomes to escape this in primordial germ cell development and early sex differenti- mechanism and invade germinal cells (Kenchington et al. ation in eukaryotes; Cg-SoxE, a homologue of Sox9,and b- 2002; Ghiselli et al. 2012). Furthermore, because ubiquitina- catenin, which are respectively expressed when sex is still not tion is also involved in the regulation of gene expression and distinguishable (or in mature females and vitellogenic oocytes plays a role in sex determination of several animals (e.g., in Crassostrea gigas [Santerre et al. 2014]); and CgFoxL2,a Caenorhabditis elegans and Drosophila;see Ghiselli et al. homologue of FoxL2, a key gene involved in ovarian determi- 2012 for a review of the literature), Ghiselli et al. (2012) pro- nation in vertebrates (Fabioux et al. 2004; Uhlenhaut and posed a relationship between ubiquitination, sex bias, and Treier 2006; Milani et al. 2011, 2015; Zhang et al. 2014). mitochondrial inheritance, and identiﬁed candidate ubiquiti- An unconventional sex-determining mechanism involving nation genes, such as the ubiquitin activating enzyme 1 (uba- doubly uniparental inheritance (DUI) of mitochondrial DNA 1) and proteasome subunit alpha 6 (psa-6), for further inves- has been hypothesized for bivalves (e.g., Breton et al. 2011; tigation. To date, this study remains the only transcriptomic Milani et al. 2016). This unique system of mitochondrial trans- analysis performed to identify genes involved in bivalve sex mission in the animal kingdom has been discovered in >100 determination and DUI. Although an in-depth study of sex- species from different bivalve orders (Mytiloida, Unionoida, speciﬁc genes has recently been conducted on male and fe- Veneroida, and Nuculanoida; Gusman et al. 2016). DUI male gonads of the gonochoric, DUI-possessing freshwater involves two distinct mitochondrial genomes that are trans- mussel Hyriopsis schlegelii (Shi et al. 2015), the link between mitted in a sex-speciﬁc way: the female-type, or F-type, trans- DUI and sex determination has not been investigated. To sum mitted through eggs and usually found in both male and up, a total of 45,422 genes differentially expressed in male female somatic tissues and in female gonadic tissues, and and female gonads in H. schlegelii and key genes reported to the male type, or M-type, transmitted through sperm and govern sex-determination pathways in mammals were iden- usually found in male gametes (reviewed in Breton et al. tiﬁed, including Sry, Dmrt1 and Sox9, upregulated in males, 2007; Passamonti and Ghiselli 2009; Zouros 2013). In addition and Foxl2 and b-catenin, upregulated in females (Shi et al. to the typical set of 13 mtDNA-encoded proteins (Gissi et al. 2015). 2008; Capt et al. 2015), novel sex-associated mtDNA- With the aim of unraveling the mechanisms underlying sex encoded proteins have been found in bivalves with DUI; F- determination and DUI in freshwater mussels, we performed 578 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE a comparative analysis of gonad transcriptome for the two sequences were also removed by specifying the appropriate DUI species Venustaconcha ellipsiformis (subfamily “ILLUMINACLIP” parameters. Reads quality was visualized Ambleminae) and Utterbackia peninsularis (subfamily with FastQC (v0.11.2; Andrews 2010). For each species, Anodontinae). Furthermore, we compared our results with high quality reads were assembled de novo to provide published sex-speciﬁc transcriptomics data from the DUI species-speciﬁc reference transcriptomes using Trinity (Haas freshwater mussel H. schlegelii and marine clams R. philippi- et al. 2013). Parameters were kept as default, except for narum and Ruditapes decussatus to identify conserved genes the –min_contig_length, which was set at 300 to remove involved in the regulation of sex-speciﬁc aspects of DUI sys- contigs shorter than 300 bp. This value was chosen for further tems. Because of the many similarities found among the annotation purposes, that is, according to the parameters distantly-related DUI species (e.g., sex-ratio bias, mitochon- ﬁxed by TransDecoder that identiﬁes ORFs with a minimum drial behavior in early embryos, rates of evolution of the length of 100 amino acids to avoid false positives (https:// two genomes; see Zouros 2013 for a review), it was indeed transdecoder.github.io/). expected that sets of genes associated with sex determination In brief, Trinity combines overlapping reads to form longer and DUI would be similar in distantly-related DUI species. We fragments called contigs, which are subsequently clustered hypothesized that if a modiﬁcation of the ubiquitination based on sequence similarity (i.e., grouping of related contigs mechanism is responsible for the retention of the paternal that correspond to portions of alternatively-spliced transcripts mtDNA in male bivalves, then “molecular signatures” of or otherwise unique portions of paralogous genes). These this modiﬁcation should be discernable in all DUI species. data were then processed to extract full-length alternatively spliced isoforms and tease apart transcripts derived from paralogous genes (Grabherr et al. 2011) to acquire nonredun- Materials and Methods dant unigenes. Sample Collection and RNA Isolation The completeness of transcriptome assemblies was tested using Benchmarking Universal Single-Copy Adult specimens of V. ellipsiformis and U. peninsularis were Orthologs (BUSCO) software (busco_v3, using the mode collected in July 2014 from Straight River (Minnesota, USA; -tran) by comparing known core eukaryotic and metazoan Lat 44.006509, Long 93.290899) and Suwannee River genes with the contigs assembled by Trinity. BUSCO pro- (Florida, USA; Lat 29.58684, Long 82.94095), respectively. vides quantitative measures for the assessment of tran- Mussels were shipped alive to the Universit e de Montr eal, scriptome completeness based on evolutionarily- where they were opened and sexed by microscopic examina- informed expectations of gene content from near- tion of gonad smears. Mature gonads of all unambiguously universal single-copy orthologs selected from different sexed specimens were frozen at 80 C until further analyses. downloadable data sets (Simao et al. 2015). Total RNA was extracted from frozen gonad tissues with RNeasy Plus Universal Mini Kit (QIAGEN Inc., Valencia, CA) and treated with Turbo DNase (AMBION, Austin, TX) follow- Annotation ing the provided protocols. A total of 25 RNA extractions were Assembled contigs were searched against two databases, performed, 14 for U. peninsularis (8 femalesand6males) and that is, Swiss-Prot and TrEMBL (Uniref90 clusters), of the 11 for V. ellipsiformis (7 females and 4 males). The quality and Universal Protein Resource (UniProt) protein database (The quantity of extracted RNA were examined via electrophoresis UniProt Consortium, 2017) withBLASTp(E-value set at on 1% agarose gel and BioDrop mLITE spectrophotometer. e05) (Altschul et al. 1990). To do this, ORFs and associ- Sixteen samples (four males and four females for each species) ated protein-coding sequences were ﬁrst predicted from were sent to the Institute for Research in Immunology and the contigs using TransDecoder (r20140704; Haas and Cancer Institute (IRIC, Universitede Montr eal, Montr eal, Papanicolaou 2017) and the longest ORF candidate was Canada) for Illumina library preparation and paired-end considered. BLASTx searches of assembled contigs against (2100 bp) sequencing (Illumina HiSeq2000; San Diego, CA). the SwissProt database were also performed to increase the proportion of annotated genes. Finally, BLASTp De Novo Transcriptome Assemblies searches for both species were performed against the Prior to assembly, raw paired-end reads were trimmed with PANMDB database. Speciﬁcally, we used the reference Trimmomatic (v0.32; Bolger et al. 2014), which is available as data set constructed by Patnaik et al. (2016), which com- part of the Trinity transcriptome assembler (trinityrna- bines protein sequence data of Arthropoda, Nematoda, seq_r20140717; Grabherr et al. 2011), to remove low- and Mollusca downloaded from the Taxonomy browser of quality reads (i.e., with Phred quality scores <30, orifthe the NCBI nr database and stored in the PANMDB (freely read contained more than ﬁve ambiguous nucleotides “N”). downloadable from the amino acid database BLAST web- As part of the trimming process, reads similar to known po- interface of the Malacological Society of Korea; http:// lymerase chain reaction primers and Illumina adapter malacol.or.kr/blast/aminoacid.html; Kang et al. 2015). Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 579 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE Read Mapping, Differential Expression, and Functional BLAST Reciprocal Analysis Enrichment Analyses A whole transcriptome reciprocal-best-BLAST-hits (RBH) anal- The de novo transcriptomes assembled for each species ysis was applied to identify orthologous genes between V. served as reference for read mapping. Speciﬁcally, raw ellipsiformis and U. peninsularis (e.g., to ﬁlter possible contam- reads were mapped using Bowtie2 (v2.2.4; Langmead ination- and species-speciﬁc transcripts). Speciﬁcally, each and Salzberg 2012), keeping the default parameters. contig from V. ellipsiformis was searched against all contigs Reads were sorted by name, and indexed with from U. peninsularis using the parameter p¼ blastnþ (default Samtools (v1.1; Li et al. 2009). Samtools idxstat was parameters were kept according to the proteinortho5.pl used to extract the numbers of reads mapped per contig script, with an E-value setate05; Lechner et al. 2011), and the results were then normalized according to and conversely each sequence of U. peninsularis was searched Robinson and Oshlack (2010). against all sequences of V. ellipsiformis. This step allowed for Two approaches were used to perform differential expres- identiﬁcation of orthologous genes between the species (e.g., sion analyses between male and female gonads for each spe- Zhao et al. 2014). From these orthologous gene sets, the cies, that is, one with the R package edgeR (version orthologs found to be differentially expressed between males edgeR_3.12.0; Robinson et al. 2010) and one with DESeq2 and females both in U. peninsularis and V. ellipsiformis using (version DESeq2_1.10.1; Love et al. 2014). We used both edgeR and DESeq2 were studied more thoroughly. approaches for comparative purposes, because multiple Transcription level and annotation were individually veriﬁed methods have been used to identify differentially expressed for each DEG in each species to conﬁrm the reliability of our genes (DEG) between male and female gonads in bivalves approach, to identify orthologous sex-biased genes that may [e.g., edgeR by Ghiselli et al. (2012) for the marine clam R. be involved in sex determination, molecular controls of DUI, or philippinarum; DESeq2 by Teaniniuraitemoana et al. (2014) both. for the oyster P. margaritifera;and Shi et al. (2015)used a statistical method developed by Audic and Claverie (1997) for Results and Discussion the freshwater mussel H. schlegelii]. DESeq2 is more sensitive Sequencing and Assemblies of the Transcriptomes to outliers and generally gives lower true positive rates, For U. peninsularis, 43 million paired end reads per individ- whereas edgeR performs better to uncover true positives ual were produced by the Illumina HiSeq2000 sequencing, (Soneson and Delorenzi 2013; Zhang et al. 2014). For a com- generating a total of 348,652,616 raw reads of 2100 bp. plete comparison of both software programs, see Zhang et al. After trimming, 288,320,738 (83%) high-quality reads were (2014). retained for transcriptome assembly. The de novo assembled Fragments per kilobase of exon per million fragment se- transcriptome using Trinity resulted in 200,961 contigs with a quenced (FPKM) were obtained using the DESeq2 N50 of 1,655 bp and average contig length of 1,026 bp Bioconductor package (Love et al. 2014), specifying (table 1). These contigs were assembled into a total of “robust¼ F” to use sum of raw counts, allowing analyses 165,788 unigenes. of differentially expressed unigenes (hereafter deﬁned as For V. ellipsiformis, 41 million reads per individual were DEG). With edgeR, DEG were selected based on FPKM values obtained on average, giving a total of 333,708,540 raw reads >1 in one or the other sex, and log2 fold change 1and generated by Illumina sequencing. After trimming, 1. Since DESeq2 removes low expression data when cal- 320,230,962 (96%) clean reads were obtained and assem- culating False Discovery Rate (FDR) values, no ﬁlter was per- bled into 285,260 contigs (N50 of 2,083 bp, average contig formed on FPKM values, and DEG were also selected based length of 1,171 bp), of which 221,362 were considered as on log2-fold change 1and 1. Only unigenes with unigenes and kept for further analyses (table 1). P-values <0.05 after adjustment by FDR (Benjamini and Assembly quality (assembled contigs) assessed by the rep- Hochberg 1995) were considered as differentially expressed. resentation of 978 core eukaryotic and metazoan genes using The package edgeR (plotSmear) was used to visualize differ- BUSCO resulted in 96.8% completely recovered genes for U. ential expression between gonad samples for each species. peninsularis (with 24.1% duplicated, 3% fragmented, and The proportion of unigenes that were differentially expressed 7% missing), and 97.8% completely recovered genes for V. according to both analytical methods was also measured. ellipsiformis (with 34.4% duplicated, 0% fragmented, and To further understand these biased expressed unigenes, 2% missing) (supplementary ﬁg. S1, Supplementary gene ontology (GO) enrichment analyses were performed Material online). In theory, the relatively high number of dupli- for each species to test whether speciﬁc GO terms (biological cates may represent alternatively spliced forms, gene duplica- process, molecular function, and cellular component levels) tion, or allelic variation (heterozygosity) in the samples used to were over-represented in DEG using GOseq (v1.28.0; construct the assemblies. Therefore, we ran the same BUSCO Young et al. 2010). An FDR corrected P-value of 0.05 was analysis using the assembled unigenes (i.e., ﬁltered for used for signiﬁcance. 580 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE Table 1 De Novo Assembly Quality Statistics of Utterbackia peninsularis and Venustaconcha ellipsiformis Transcriptomes U. peninsularis V. ellipsiformis Total raws 348,652,616 333,708,540 Average raws 43,581,577 41,713,567 Total trimmed raws 288,320,738 320,230,962 Average trimmed raws 36,040,092 40,028,870 Total trinity “unigenes” 165,788 221,362 Total trinity contig 200,961 285,260 Percent GC 33.36 37.58 Contig N10 5,581 6,494 Contig N20 3,895 4,650 Contig N30 2,916 3,560 Contig N40 2,224 2,731 Contig N50 1,655 2,083 Median contig length 558 595 Average contig 1,026 1,171 Total assembled bases 206,277,984 334,013,631 Unigene N10 4,882 5,792 Unigene N20 3,298 3,896 Unigene N30 2,409 2,789 Unigene N40 1,750 2,019 Unigene N50 1,251 1,433 Median unigene length 510 516 Average unigene 884 944 Total assembled bases 146,628,566 208,707,657 alternatively spliced forms and paralogs) to better understand the origin of our duplicates (supplementary ﬁg. S1, Supplementary Material online). This analysis still showed a similar number of completely recovered genes (96.6% for U. peninsularis and 97.5% for V. ellipsiformis), but a markedly lower number of duplicates (13.9% and 24.3%, respectively), suggesting that gene duplication and/or alternative splicing may contribute to the relatively high number of duplicates. Nevertheless, the high number of complete genes that was recovered provides an important validation of the complete- ness of both assemblies. FIG.1.—Taxonomic distribution of top BLASTp hits in Uniref90 pro- For both species, the raw sequence data in FASTQ format tein database for U. peninsularis contigs (A)and V. ellipsiformis contigs (B). have been submitted to the National Centre for Numbers of unique hits in each group are shown. Biotechnology Information (NCBI) Sequence Read Archive (SRA) database and are accessible under accession numbers SRR6279376-83 (U. peninsularis) and SRR6279374-75 and taxon Mammalia (supplementary ﬁg. S2, Supplementary SRR6279384-89 (V. ellipsiformis). Material online). With Uniref90, the largest proportion of hits (46.4%) came from the oyster species C. gigas Sequence Annotation (ﬁg. 1), for which whole-genome sequencing and anno- BLASTp searches against SwissProt and Uniref90, BLASTx tation were recently completed (Zhang et al. 2012). The against SwissProt, and BLASTx against PANMDB, respec- second largest proportion of hits (23.6%) came from the tively, allowed annotating 26,469 (13.1%), 32,515 molluscan class Gastropoda (ﬁg. 1). These results are com- (16.2%), 28,495 (14.2%), and 32,503 (16.1%) contigs parable with those reported in other recent de novo tran- for U. peninsularis. All databases together gave a total scriptome sequencing studies for freshwater mussels or annotation success of 16.8% contigs. With SwissProt, other bivalve species (i.e., 15–22% annotated contigs/ most of the matches came from the extensively studied unigenes; Shi et al. 2015; Patnaik et al. 2016; Yarra Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 581 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE et al. 2016; Wang et al. 2017) and indicate that it can be total, we retained 37,947 genes with signiﬁcant blast hits challenging to provide good annotations for nonmodel in our downstream analyses. organisms for which few genomic resources are available. Also, very few contigs matched bacterial sequences (or DEG between Male and Female Gonads in Each Species other possible contaminating species such as trematodes), suggesting negligible contamination for U. peninsularis, The differential gene expression analysis between sexes for U. and indicating that the assembly for this species com- peninsularis revealed 7,281 DEG with edgeR, and 23,107 prised mainly gonadal coding RNAs. DEG with DESeq2. The DEG identiﬁed using DESeq2 con- These results obtained for U. peninsularis, however, are in tained all the DEG identiﬁed with edgeR (ﬁg. 2). A total of sharp contrast with those obtained for V. ellipsiformis,for 4,315 unigenes were shown to be upregulated in males which 63,126 (22.1%) contigs were annotated with versus 2,966 in females with edgeR, whereas 14,216 genes Uniref90 but with almost half of them matching trematode were considered upregulated in males versus 8,891 in females species (ﬁg. 1). The fractions of contigs annotated with with DESeq2. This is a smaller number than what was found SwissProt and PANMDB were also slightly larger for V. ellipsi- in the DUI-possessing freshwater mussel H. schlegelii using a formis, with 49,904 (17.5%) and 50,444 (17.7%) contigs statistical text proposed by Audic and Claverie (1997) (45,422 annotated with BLASTp and BLASTx searches against DEG, with 19,511 unigenes upregulated in male gonads vs. SwissProt (see supplementary ﬁg. S2, Supplementary 25,911 in female gonads; Shi et al. 2015), but a substantial Material online), and 57,215 (20%) contigs annotated with number considering what has been revealed for example in PANMDB. Again, several matches came from C. gigas with the oyster species Pinctada margaritifera using DESeq2 (1,993 Uniref90 (ﬁg. 1), but the discovery that 40% of the anno- DEG, with 1,419 contigs upregulated in male gonads vs. 574 tated contigs had a signiﬁcant match with trematode species, in female gonads; Teaniniuraitemoana et al. 2014)and in the which is not thecasefor U. peninsularis, strongly suggest a DUI-possessing marine clam R. philippinarum using edgeR trematode infection in V. ellipsiformis, a phenomenon that (1,575 DEG, with most of the genes being male biased; is often observed in freshwater mussels, which are inter- Ghiselli et al. 2012). Although different DEG analytical meth- mediate hosts of several trematode species (Gangloff et al. ods and experimental setups may explain these dissimilarities, 2008). This contamination most likely explains the higher it was suggested that a low number of DEG should be number of duplicates obtained with BUSCO as well as the expected in bivalves considering their general lack of second- higher percentage of annotated contigs for V. ellipsifor- ary sexual characters and sexual dimorphism, and also be- mis. We performed additional analyses to verify the extent cause sex-speciﬁc function of reproductive genes in these of this contamination in our samples, and found that three organisms is limited to gonad development, gametogenesis, of four female samples, but no male samples, produced and fertilization (Ghiselli et al. 2012). However, contrary to several hits when blasted against the three most expressed most bivalve species, unionoid freshwater mussels have a trematode transcripts (data not shown). This is consistent unique reproductive mode that may cause pronounced sexual with recent report that mature female freshwater mussels dimorphism: 1) their larvae mature in certain areas unique to are more susceptible to trematode infection (Mu¨ lleretal. female gills called marsupia (which can also alter shell mor- 2015). phology) and 2) because most species require a host ﬁsh for We discarded contaminant trematode sequences from their larvae, many gravid females, including those of V. ellipsi- the transcriptome of V. ellipsiformis by relying on annota- formis (Allen et al. 2007), display adaptations such as mantle- tion to remove Platyhelminthes transcripts (e.g., Sayadi derived lures to actively attract a prospective ﬁsh for their par- et al. 2016). Although this approach may result in the sig- asitic larvae (Haag and Warren 1999; Barnhart et al. 2008; niﬁcant loss of mussel-speciﬁc data, we considered it to be a Zieritz and Aldridge 2011; Haag 2012). Interestingly, sexual valid approach for retaining the maximum number of dimorphism in burrowing behavior has also been reported in sequences related to our focal species (V. ellipsiformis). the freshwater mussel DUI species Elliptio complanata (Flynn Alternatively, this approach can lead to poor decontamina- et al. 2013). To our knowledge, however, no sexual dimor- tion, that is, with many nonannotated contigs belonging to phism has been reported in U. peninsularis. Mostofthe DEG nontarget or contaminant species. However, because the found in this species are male-biased, which is consistent with ultimate goal of this study was to identify similar sets of what has been found in most bivalve species, as well as in many genes associated with sex determination and DUI in two other animals (Ghiselli et al. 2012; Teaniniuraitemoana et al. different freshwater mussel species, we tried to minimize 2014; Peng et al. 2015; Li et al. 2016), and may be explained by the contamination problem by using the best matching re- a higher proportion of essential functions for female-biased ciprocal BLAST between both species to ﬁlter out genes (those shared by the two sexes) than for male-biased contamination-speciﬁc (and species-speciﬁc) contigs (see genes (see Ghiselli et al. 2012). This explanation is also consis- below). A total of 25,179 (40%) annotated contigs were tent with the previous assumption that male development is considered as probable contaminants and discarded. In associated with an activation of several testis-speciﬁc genes 582 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE FIG.2.—Smear plots of DEG between male and female gonads in U. peninsularis (A)and V. ellipsiformis (B). Black dots represent genes normally expressed in both tissues. Green dots, DEG revealed by edgeR software; orange dot, DEG revealed by DESeq2. Green dots surrounded by orange represent genes considered as DEG by both software packages. Positive fold-changes represent genes upregulated in males, while negative ones represent genes upregulated in females. and/or a repression of several genes vital for ovarian develop- revealed by edgeR (765 upregulated in males and 9,004 in ment (see Peng et al. 2015). females) and 52,652 revealed by DESeq2 (3,317 upregulated For V. ellipsiformis, 11,408 DEG were revealed with edgeR in males and 49,335 in females). These numbers, although and 52,257 DEG with DESeq2 (ﬁg. 2). Only 60 unigenes were not surprising since the percentage of annotation (< 22%), upregulated in males versus 11,348 in females with edgeR, was very low, either suggest 1) that the contamination could whereas DESeq2 revealed 3,134 unigenes upregulated in not be completely removed from our data using the approach males and 49,123 in females. This important female-biased described above, or 2) that changes in female host gene ex- gene expression in V. ellipsiformis is most likely explained by pression occurred in response to parasite infection (e.g., trematode contamination found in three out of four female Wijayawardena et al. 2016). Fortunately, because our U. pen- gonad samples (see the annotation section above). After hav- insularis samples were not infected by trematodes, our BLAST ing ﬁltered out 25,179 Platyhelminthes contigs from our data, reciprocal analysis discussed below allowed us to mitigate this the number of DEG did not vary much with 9,769 DEG problem and to accomplish our main objective, that is, identify Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 583 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE FIG.3.—Signiﬁcantly enriched-GO terms involved in molecular function (top), biological process (middle), and cellular component (bottom) in ovary (pink) and testis (blue) of both species. Analyses were performed on DEGs from blast reciprocal analyses between V. ellipsiformis and U. peninsularis. X-axis shows the number of DEGs contained in each GO category. similar sets of genes associated with sex determination and signiﬁcantly enriched in both species (supplementary table S2, DUI in two different freshwater mussel species. Supplementary Material online). Top GO terms enriched in female and male gonads of both species are shown in ﬁgure 3 (for complete GOSeq results for U. peninsularis and V. ellipsi- BLAST Reciprocal Analysis and Enriched GO Terms and formis after trimming for trematode contigs, see supplemen- Pathways in Male and Female Gonads tary tables S3 and S4, respectively, Supplementary Material online). In transcriptomes of both species, GO terms show To identify the main genes possibly involved in sex determi- higher counts for male-biased genes. At the biological pro- nation and/or DUI in U. peninsularis and V. ellipsiformis,we cess level, single-organism cellular process, cellular compo- performed a whole transcriptome BLAST analysis. From the nent organization or biogenesis, and organelle organization reciprocal BLAST results, we selected the genes found to be had the top counts. The intracellular and organelle components differentially expressed in U. peninsularis in our previous anal- represented the majority of terms at the cellular component yses using edgeR and DESeq2, and veriﬁed their sex-bias in V. level, whereas at the molecular function level, protein binding ellipsiformis DEG ﬁltered for “ﬂatworm contigs” (i.e., we ver- was the most represented term (ﬁg. 3). Similar results were also iﬁed that every male-biased or female-biased gene in one obtained for other bivalve species, including the oyster species was also male-biased of female-biased in the other). Crassostrea hongkongensis (Tong et al. 2015), the scallop A total of 12,000 orthologous relationships was found Patinopecten yessoensis (Li et al. 2016), as well as the DUI- between both species, 16,740 when isoforms were included. containing marine clam R. philippinarum (Ghiselli et al. 2012) From these orthologous genes, 2,342 were found to be dif- and freshwater mussel Cristaria plicata (Patnaik et al. 2016). ferentially expressed in U. peninsularis by edgeR and 2,583 Unexpectedly, “reproduction” and “ubiquitination,” the two (including the 2,342 from edgeR) by DESeq2. Speciﬁc atten- categories on which Ghiselli et al. (2012) previously focused be- tion was given to these 2,583 genes (1,568 male-biased and cause of their direct relationship with sex determination and DUI, 1,015 female-biased) to increase the probability of identifying were not signiﬁcantly overrepresented in our results. We thus candidate genes that could have a role in sex determination concentrated our effort on the reciprocal BLAST/DEG analyses. and/or DUI. Of the 1,568 male-biased orthologous genes, 1,432 (91.3%) had a signiﬁcant BLASTp (or BLASTx) match in the SwissProt and/or UniRef90 protein databases (E-value Identiﬁcation of Candidate Genes Involved in Sex set at e05), whereas of the 1,015 female-biased ortholo- Determination gous genes, 903 (89%) had a signiﬁcant BLASTp (or BLASTx) match in the SwissProt and/or UniRef90 protein databases To identify candidate genes that contribute to sex deter- (supplementary table S1, Supplementary Material online). mination/differentiation, or maintenance in mature The same genes were also found to be differentially expressed gonads, we ﬁrst searched in our reciprocal BLAST/DEG in V. ellipsiformis. results for sex-speciﬁc expression patterns (i.e., tran- Using all annotated and differentially expressed ortholo- scribedonly inone sex)—asexpectedfor sex- gous genes, GOSeq analyses revealed 57 ontology categories determining genes (table 2). In addition, we identiﬁed 584 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE Table 2 Sex-Speciﬁcally Expressed Genes in Male and Female Gonads of Utterbackia peninsularis and Venustaconcha ellipsiformis Unigene ID Gene Symbol Description Sex (fold change) c43183_g1 – Uncharacterized protein M (10) c63770_g1 – Uncharacterized protein M (10) c106662_g8 TSSK3_MOUSE Testis-speciﬁc serine/threonine-protein kinase 3 M (11) c69491_g1 RDS2_CHICK Photoreceptor outer segment membrane glycoprotein 2 M (12) c106233_g1 EDR1_ARATH Serine/threonine-protein kinase EDR1 M (12) c55198_g1 CE049_MOUSE Uncharacterized protein C5orf49 homolog M (12) c85464_g1 CALL3_MOUSE Calmodulin-like protein 3 F (9) c66486_g1 SMCO3_MOUSE Single-pass membrane and coiled-coil domain-containing protein 3 F (13) c237103_g1 – Uncharacterized protein F (13) Note.—Unigene ID indicates the names of the unigene given by Trinity to the assembly ﬁles. Gene symbol indicates the annotation format provided by the Uniprot database followed by a quick short description. Sex indicates in which sex the unigene is speciﬁcally expressed, fold change values (obtained from DESEq2 software) are based on the results obtained for the noncontaminated U. peninsularis. Table 3 Top Ten Genes Showing the Greatest Difference in Expression in Male and Female Gonads of Utterbackia peninsularis and Venustaconcha ellipsiformis Unigene ID Gene Symbol Description Sex (fold change) c1056_g1 K1S6Q8_CRAGI Serine/threonine-protein phosphatase 6 regulatory ankyrin repeat subunit B M (12) c235886_g1 KLH10_HUMAN Kelch-like protein 10 M (11) c104796_g5 K1RU43_CRAGI WD repeat-containing protein on Y chromosome M (11) c47982_g1 SAM15_MACFA Sterile alpha motif domain-containing protein 15 M (11) c847_g1 K1Q7E2_CRAGI Sperm motility kinase X M (11) c46585_g1 K1RFU6_CRAGI Proteasome activator complex subunit 3 M (11) c173681_g1 K1RFD4_CRAGI Probable 4-coumarate–CoA ligase 3 M (11) c44943_g1 PYG_DROME Glycogen phosphorylase M (11) c36459_g1 CC151_BOVIN Coiled-coil domain-containing protein 151 M (11) c110789_g1 K1R8H1_CRAGI Testis-speciﬁc serine/threonine-protein kinase 5 M (11) c69487_g1 LECG_THANI Galactose-speciﬁc lectin nattectin; F (14) c92390_g1 LECG_THANI Galactose-speciﬁc lectin nattectin; F (14) c82290_g1 LECM2_ERYPO C-type lectin lectoxin-Lio2; F (14) c83693_g1 PSM_MYTCA Shell matrix protein; F (14) c79356_g1 LECM1_PHIOL C-type lectin lectoxin-Phi1; F (14) c99224_g1 PLC_HALLA Perlucin; F(13) c173670_g1 LPSBP_PERAM Hemolymph lipopolysaccharide-binding protein; F (13) c96153_g1 LECG_THANI Galactose-speciﬁc lectin nattectin; F (13) c105339_g1 LECM2_ERYPO C-type lectin lectoxin-Lio2; F (13) c55454_g1 CL17A_HUMAN C-type lectin domain family 17, member A; F (13) Note.—Unigene ID indicates the names of the unigene given by Trinity to the assembly ﬁles. Gene symbol indicates the annotation format provided by the SwissProt or Uniprot database followed by a short description. Sex indicates in which sex the unigene is upregulated, fold change values (obtained from DESEq2 software) are based on the results obtained for the noncontaminated U. peninsularis the top ten upregulated (and annotated) genes in male the ubiquitination and subsequent proteasomal degrada- and female gonads of both freshwater mussel species (ta- tion of target proteins during spermatogenesis (Arama ble 3). Among these, six genes were exclusively expressed et al. 2007). In addition, sperm motility kinase 2B and in male gonads and three in female gonads. Not surpris- WD repeat-containing protein on Y chromosome are ingly, these sex-speciﬁcally expressed and upregulated involved in sperm fertility and motility, whereas testis- genes appear to play roles in spermatogenesis and oogen- speciﬁc serine/threonine-protein kinases 3 and 5 are in- esis, which are the ultimate biological processes in the volved in germ cell survival during meiosis (e.g., Peng male and female gonads, and in reproduction in general. et al. 2015). Among other male-biased sequences, we In males, for example Kelch-like protein 10 and also identiﬁed genes homologous to spermatogenesis- Proteasome activator complex are members of the E3 associated protein (SPATA 4, 17, 22, and 24), sperm ﬂa- ubiquitin–protein ligase complex process, which mediates gellar protein 1/2, sperm motility kinase, sperm-associated Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 585 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE antigen and surface protein,and spermine oxidase poten- Beta-catenin homolog is named Armadillo andisknownto tially expressed in the male reproductive tissues (supple- be involved (among other processes) in oogenesis and ovarian mentary table S4, Supplementary Material online). follicle cell development. The Rspo1/Wnt/Beta-catenin path- Speciﬁcally expressed and top upregulated genes in way isalsoinvolvedinbird sex determination(Beukeboom females were mostly lectin-type genes, which are often highly and Perrin 2014). represented in ovary transcriptomes in animal species, includ- Collectively, the results presented in table 4 suggest shared ing bivalves. These genes have a number of biological func- sex-determining pathway genes in distantly-related, as well as tions such as immune defense, block to polyspermy at in SMI and DUI bivalves. Zhang et al. (2014) suggested a fertilization, and sperm–egg recognition (Luckenbach et al. deeply conserved role in sex determination for Dmrt1/Dsx/ 2008; Moy et al. 2008; Patnaik et al. 2016). Among other MAB-3 in both invertebrates and vertebrates, but also for female-biased genes, we identiﬁed vitellogenin, putative vitel- Sry and FoxL2, even if they were thought to be new recruits logenin receptor,and Hsp90 organizing protein,which are to sex-determining pathways in vertebrates or placental mam- involved in vitellogenesis, a central process in oogenesis. We mals (Gamble and Zarkower 2012; Matson and Zarkower also found cyclin-dependent kinase 2-associated protein 2, 2012). Indeed, some steps of the sex determination pathway which is involved in another vital process in oogenesis, that may be deeply conserved in metazoans, despite rapid evolu- is, oocyte maturation, and superoxide dismutase [Cu–Zn], tion of the regulatory pathways that in bivalves, for example, which could neutralize reactive oxygen species and ensure may involve both genetic and environmental factors (Zhang oocyte quality (Peng et al. 2015 and above). et al. 2014; Breton et al. 2017). To better identify genes related to sex-determining path- Zhang et al. (2014) proposed a working model for sex de- ways, we also searched in our reciprocal BLAST/DEG results termination in the oyster C. gigas,in which CgSoxH,a Sry-like for candidate genes previously identiﬁed in the literature in 1) gene that is strictly expressed in testis, would play a leading the model organisms Ca. elegans, Drosophila melanogaster, role in the sex-determining pathway by directly or indirectly Mus musculus (see Zhang et al. 2014 for the list); 2) four activating CgDsx, a DM domain gene like those (e.g., Dmrt1) bivalve species without DUI, that is, the oysters C. gigas and that have been identiﬁed as master switches for testis devel- C. honkongensis, the scallop P. yessoensis and the clam R. opment in animal species studied so far. Both CgSoxH and decussatus (Zhang et al. 2014; Tong et al. 2015; Li et al. CgDsx would interact with, or inhibit, CgFoxL2, which is usu- 2016; Ghiselli F, Iannello M, Puccio G, Chang PL, Plazzi F, ally speciﬁcally expressed in ovaries (this scenario is consistent Nuzhdin SV, Passamonti M. unpublished data); and 3) two with the reported interaction among Sry, Dmrt1,and FoxL2 in bivalve species with DUI, that is, the marine clams R. philip- mammals). Even if the number of species studied is still lim- pinarum and the freshwater mussel H. schlegelii (Ghiselli et al. ited, we suggest that this general model probably applies to 2012; Shi et al. 2015). Of the 26 key sex-determining pathway most bivalve species. We also propose that major interspeciﬁc genes examined, putative homologs were found for 8 genes differences in bivalves will be found in cis/trans genetic ele- or gene families in our two freshwater mussel species (table 4; ments and/or environmental factors that may control these supplementary table S5, Supplementary Material online). Of three “sex-determining key genes.” For example, a mito- these, four genes (Sry, Fem-1, MAB-3,and FoxL2)were dif- nuclear sex determination system in which mitochondrially ferentially expressed but did not show sex-speciﬁc expression encoded elements distorting sex ratios has been hypothesized in our samples, although we cannot rule out the possibility to occur in species with DUI (e.g., Breton et al. 2011; Milani that they may have sex-speciﬁc expression at an earlier devel- et al. 2016; Pozzi et al. 2017). In both plants and animals, opmental stage. For example, Sry/Sox30 was found to be cases of mitochondrial genomes acting on germ line develop- expressed only in male oysters and clams (table 4). ment and/or with sex ratio distortion properties have been As previously mentioned, both Sry (sex determining region- documented (e.g., Chase 2007; Perlman et al. 2015). In y)-box 30 (Sox 30)and Dmrt1 (MAB-3 and dsx related tran- bivalves with DUI, the existence of two sex-speciﬁc mitochon- scription factor 1) are known to be involved in male sex drial genomes on which selection can operate represents a determination in nematodes, fruit ﬂies, and vertebrates possible source of novel sex determination mechanisms. (Wallis et al. 2008; Kopp 2012), whereas FoxL2 is known for its role in ovarian determination in vertebrates (Uhlenhaut and Identiﬁcation of Candidate Genes Involved in DUI Treier 2006). Fem-1 encodes an ankyrin repeat-containing pro- tein orthologous to human FEM1A and is involved in male sex- As mentioned earlier, a modiﬁcation of the ubiquitin–protea- determination in Ca. elegans, and also in the protein ubiquiti- some system has been hypothesized as responsible for the nation pathway (Doniach and Hodgkin 1984). Two other retention of sperm mitochondria in male embryos of species genes that seem well represented in bivalves, and more gen- with DUI (Kenchington et al. 2002; Ghiselli et al. 2012). To erally in mollusks (Kim et al. 2017), are Bar-1/Arm/Ctnnb1 identify putative factors that could be involved in the degrada- (Beta-catenin)and Tra-1 (table 4). Tra-1 controls female so- tion of sperm mitochondria through ubiquitination in DUI matic sex differentiation in Drosophila,whereas Drosophila bivalves, Milani et al. (2013) analyzed structure and localization 586 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 587 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Table 4 Key Sex-Determining Pathway Genes in Model Organisms and Bivalve Species Gene Caenorhabditis Drosophila Mus Crassostrea Crassostrea Patinopecten Ruditapes Ruditapes Hyriopsis Utterbackia elegans melanogaster musculus gigas hongkongensis yessoensis decussatus philippinarum schlegelii peninsularis/ Venustaconcha ellipsiformis Amh (Mmu) ms Amhr2 (Mmu), Bar-1 (Cel), Arm (Dme), Ctnnb1 (Mmu) fs fb Cbx2 (Mmu) ms Doa (Dme) fs Dsx (Dme), Dmrt1 (Mmu), Mab-3 (Cel) ms ms ms ms ms mb mb mb Fog-1,-2,-3 (Cel, Dme, Mmu) ms ms Fem-1 (Cel, Dme), Fem1b (Mmu) ms mb mb mb FoxL2 (Mmu) fs fs fs fb fb fb Fru (Dme) ms Zglp1 (Mmu) Her-1 (Cel), Her (Dme) fs Nr0b1 (Mmu) ms Rspo1 (Mmu) fs fb Rnt-1 (Cel), Run (Dme) fs Sdc (Cel, Dme, Mmu) fs Sxl (Dme) fs Sry/Sox30/SoxH (Mmu) ms ms ms mb ms ms mb mb Sis-a (Dme) fs Sf1 (Mmu) ms Sox9 (Mmu), Sox100B (Dme) ms mb Gata4 (Mmu) mb mb Tra-1-2 -3(Cel), Tra2 (Dme) fs fs fb Wt1 (Mmu) ms mb Wnt4 (Dme, Mmu) fs fb Xol-1 (Cel) ms Note.—Names in parentheses indicate the species in which each gene has been identiﬁed; shaded box indicates presence of the gene in the genome/transcriptome; fs, female-speciﬁc expression (i.e., transcribed only in females); ms, male-speciﬁc expression; fb, female-biased expression; mb, male-biased expression. See supplementary table S5, Supplementary Material online, for a short description of each gene Capt et al. GBE of three transcripts—that is, baculoviral IAP repeat-containing Grabowy, 1983; Umen and Goodenough 2001). A growing 4(birc), proteasome subunit alpha 6 (psa), and AN1 zinc ﬁnger body of literature has also shown that epigenetic modiﬁca- ubiquitin-like domain (anubl1)—previously shown to present tions, such as DNA methylation, can also regulate sex deter- sex-speciﬁc and family biases (i.e., family producing predom- mination and differentiation (see Liu et al. 2015). Strikingly, in inantly females or males) in the Manila clam R. philippinarum. the leafhopper Zyginidia pullula, the symbiont Wolbachia In situ hybridization conﬁrmed the localization of these pipientis disrupts male imprinting by modifying methylation transcripts in gametogenic cells (Milani et al. 2013). Also, patterns in gonads; as a result, Wolbachia feminizes genetic homologs of these genes were shown to be involved in males of Z. pullula (Negri et al. 2009). In our results, transcripts reproduction and ubiquitination in other animals: PSA is a associated with DNA methyltransferases, histone deacety- subunit of the proteasome, a complex known to be involved lases, and histone acetyltransferases were identiﬁed, and in male sexual differentiation, whereas ANUBL1 possesses an most showed male-biased expression (supplementary table ubiquitin domain, and BIRC acts as an inhibitor of apoptosis S1, Supplementary Material online). Among these male- (Milani et al. 2013). For these reasons, the authors hypothe- biased transcripts, DNMT1, a DNA methyltransferase associ- sized that these genes could have a role in sex determination ated with mitochondria, was detected. However, because the and could also be responsible for the maintenance/degrada- epigenetic mechanisms underlying organellar DNA inheri- tion of spermatozoon mitochondria during embryo develop- tance and/or sex determination/differentiation are still poorly ment of the DUI species R. philippinarum.It is worth understood, we cannot infer any detailed functions of these mentioning that PSA hasalsobeenreportedinmale-biased genes from their expression patterns. Future studies, currently families in the DUI species Mytilus edulis (Diz et al. 2013), underway in our laboratory, are clearly needed to reveal the supporting the hypothesis of an involvement of the ubiquiti- complex interplay between methylation, sex bias, and mito- nation system (and of the proteasome subunit alpha) in sex chondrial inheritance in bivalves with DUI. determination in bivalves. To identify genes related to the DUI mechanism, we Conclusion searched our reciprocal BLAST results for genes previously In this study, we present the ﬁrst comparative gonadal tran- suggested to be involved in this process in the DUI species scriptomic analysis of two freshwater mussel species aimed at R. philippinarum and M. edulis, that is, for genes involved in assessing the link between DUI of mtDNA and sex determi- ubiquitination that show male-biased transcription, indicating nation in Unionoida. Using a BLAST reciprocal analysis to iden- their potential role in spermatogenesis (Ghiselli et al. 2012; Diz tify orthologs between V. ellipsiformis and U. peninsularis and et al. 2013; Milani et al. 2013). Annotation, function, and comparing our results with previously published sex-speciﬁc expression ratio of male-biased (and female-biased) ubiquiti- transcriptomics data from distantly-related DUI species, we nation genes are reported in supplementary table S1, were able to identify elements possibly involved in the regu- Supplementary Material online. Compared to females in which lation of sex-speciﬁc aspects of DUI systems. Our search for 18 ubiquitination genes were differentially expressed, 53 ubiq- candidate genes involved in sex determination suggests that uitination genes were male-biased. Among them, psa6 and some steps of the sex determination pathway may be deeply anubl1, two male-biased ubiquitination genes were absent conserved in metazoans, whereas our search for candidate from female-biased genes, as were birc5 and birc2.These genes involved in DUI supports previous hypotheses that a results, consistent with the hypothesis of Milani et al. (2013), modiﬁcation of the ubiquitination mechanism could be re- suggest thatpsa6 might be important in sex determination and/ sponsible for the retention of the paternal mtDNA in male or differentiation in DUI species, that anubl1 could be a factor bivalves. Our results also suggest that DNA methylation could tagging sperm mitochondria that differentiate them from egg be involved in the maintenance of DUI in bivalves. mitochondria, and that birc could be an additional mitochon- drial tag that would protect the midpiece mitochondria from degradation. Our results are also consistent with our hypothesis Supplementary Material that sets of genes associated with sex determination and DUI Supplementary data are available at Genome Biology and are similar in distantly-related DUI species, as exempliﬁed here Evolution online. with freshwater mussels versus a marine clam and a marine mussel. However, additional RNA- and protein-based studies such as those conducted by Milani et al. (2013) and Diz et al. Acknowledgments (2013) will be needed to conﬁrm our results. Lastly, although speculative, one of the factors that This study was supported by Natural Sciences and Engineering could also inﬂuence inheritance of mitochondrial genes Research Council Discovery Grants awarded to S.B. (RGPIN/ is DNA methylation. For example, in the green alga 435656-2013) and D.T.S. (RGPIN/217175-2013). Any use of Chlamydomonas, maternal inheritance of chloroplast DNA is trade, ﬁrm, or product names is for descriptive purposes only thought to be inﬂuenced by DNA methylation (Sager and and does not imply endorsement by the US Government. 588 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Link between DUI of mtDNA and Sex Determination in Bivalves GBE bias, mitochondrial doubly uniparental inheritance and sex determina- Literature Cited tion. Mol Biol Evol. 29(2):771–786. Allen DC, et al. 2007. Early life-history and conservation status of Gissi C, Iannelli F, Pesole G. 2008. Evolution of the mitochondrial genome Venustaconcha ellipsiformis (Bivalvia, Unionidae) in Minnesota. Am of Metazoa as exempliﬁed by comparison of congeneric species. Midl Nat. 157:74–91. Heredity 101(4):301–320. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local Grabherr MG, et al. 2011. Full-length transcriptome assembly from RNA- alignment search tool. J Mol Biol. 215(3):403–410. Seq data without a reference genome. Nat Biotechnol. Andrews S. 2010. FastQC: a quality control tool for high throughput se- 29(7):644–652. quence data. Babraham Bioinform. Available from: http://www.bioin- Gusman A, Lecomte S, Stewart DT, Passamonti M, Breton S. 2016. Pursuing formatics.babraham.ac.uk/projects/fastqc/ the quest for better understanding the taxonomic distribution of the Arama E, Bader M, Rieckhof GE, Steller H. 2007. A ubiquitin ligase com- system of doubly uniparental inheritance of mtDNA. PeerJ 4:e2760. plex regulates caspase activation during sperm differentiation in Haag WR. 2012. North American freshwater mussels: natural history, ecol- Drosophila. PLoS Biol. 5(10):e251. ogy, and conservation.Host use and host infection strategies. Audic S, Claverie J-M. 1997. The signiﬁcance of digital gene expression Cambridge University Press. p. 140–179. proﬁles. Genome Res. 7(10):986–995. Haag WR, Warren ML. 1999. Mantle displays of freshwater mussels elicit Barnhart CM, Haag WR, Roston WN. 2008. Adaptations to host infection attacks from ﬁsh. Fresh Biol. 42(1):35–40. and larval parasitism in Unionoida. J N Am Benthol Soc. Haas BJ, Papanicolaou A. 2017. TransDecoder. Retrieved from http://trans- 27(2):370–374. decoder.github.io Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a Kang SW, et al. 2015. Construction of PANM database (Protostome DB) practical and powerful approach to multiple testing. J R Stat Soc for rapid annotation of NGS data in Mollusks. Korean J Malacol. Series B Stat Methodol. 57:289–300. 31(3):243–247. Beukeboom LW, Perrin N. 2014. The evolution of sex determination. Kenchington E, MacDonald B, Cao L, Tsagkarakis D, Zouros E. 2002. Oxford: University Press. Genetics of mother-dependent sex ratio in blue mussels (Mytilus Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a ﬂexible trimmer for spp.) and implications for doubly uniparental inheritance of mitochon- Illumina sequence data. Bioinformatics 30(15):2114–2120. drial DNA. Genetics 161(4):1579–1588. Breton S, et al. 2009. Comparative mitochondrial genomics of freshwater Kim MA, et al. 2017. Alternative splicing proﬁle and sex-preferential gene mussels (Bivalvia: Unionoida) with doubly uniparental inheritance of expression in the female and male paciﬁc abalone Haliotis discus han- mtDNA: gender-speciﬁc open reading frames and putative origins of nai. Genes (Basel) 8(3):99. replication. Genetics 183(4):1575–1589. Kopp A. 2012. Dmrt genes in the development and evolution of sexual Breton S, et al. 2011. Novel protein genes in animal mtDNA: a new sex dimorphism. Trends Genet. 28(4):175–184. determination system in freshwater mussels (Bivalvia: Unionoida)? Mol Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie Biol Evol. 28(5):1645–1659. 2. Nat Rev Clin Oncol. 9(4):357–359. Breton S, Beaupre HD, Stewart DT, Hoeh WR, Blier PU. 2007. The unusual Lechner M, et al. 2011. Proteinortho: detection of (Co-)orthologs in large- system of doubly uniparental inheritance of mtDNA: isn’t one scale analysis. BMC Bioinformatics 12:124. enough? Trends Genet. 23(9):465–474. Li H, et al. 2009. The sequence alignment/map format and SAMtools. Breton S, Capt C, Guerra D, Stewart D. 2017. Sex determining mecha- Bioinformatics 25(16):2078–2079. nisms in bivalves. Preprints 2017:2017060127. Li Y, et al. 2016. Transcriptome sequencing and comparative analysis of Capt C, Passamonti M, Breton S. 2015. The human mitochondrial genome ovary and testis identiﬁes potential key sex-related genes and path- may code for more than 13 proteins. Mitochondrial DNA ways in scallop Patinopecten yessoensis. Mar Biotechnol. 27:3098–3101. 18(4):453–465. Chase CD. 2007. Cytoplasmic male sterility: a window to the world of Liu H, et al. 2015. Large-scale transcriptome sequencing reveals novel plant mitochondrial–nuclear interactions. Trends Genet. 23(2):81–90. expression patterns for key sex-related genes in a sex-changing ﬁsh. Chavez-Villalba J, et al. 2011. Determination of gender in the pearl oyster Biol Sex Differ. 6(1):26. Pinctada margaritifera. J Shellﬁsh Res. 30(2):231–240. Llera-Herrera R, Garcıa-Gasca A, Abreu-Goodger C, Huvet A, Ibarra AM. Coe WR. 1943. Sexual differentiation in Mollusks. I. Pelecypods. Q Rev 2013. Identiﬁcation of male gametogenesis expressed genes from the Biol. 18(2):154–164. scallop Nodipecten subnodosus by suppressive subtraction hybridiza- Diz AP, et al. 2013. Proteomic analysis of eggs from Mytilus edulis females tion and pyrosequencing. PLoS One 8(9):e73176. differing in mitochondrial DNA transmission mode. Mol Cell Love MI, Huber W, Anders S. 2014. Moderated estimation of fold Proteomics 12(11):3068–3080. change and dispersion for RNA-seq data with DESeq2. Genome Biol. Doniach T, Hodgkin J. 1984. A sex-determining gene, fem-1,required for 15(12):550. both male and hermaphrodite development in Caenorhabditis ele- Luckenbach JA, Iliev DB, Goetz FW, Swanson P. 2008. Identiﬁcation of gans. Dev Biol. 106(1):223–235. differentially expressed ovarian genes during primary and early sec- Fabioux C, et al. 2004. Oyster vasa-like gene as a marker of the germline ondary oocyte growth in coho salmon, Oncorhynchus kisutch.Reprod cell development in Crassostrea gigas. Biochem Biophys Res Commun. Biol Endocrinol. 6:2. 320(2):592–598. Matson CK, Zarkower D. 2012. Sex and the singular DM domain: insights Flynn K, et al. 2013. Burrowing in the freshwater mussel Elliptio compla- into sexual regulation, evolution and plasticity. Nat Rev Genet. nata is sexually dimorphic and feminized by low levels of atrazine. J 13(3):163–174. Toxicol Environ Health A 76(20):1168–1181. Milani L, Ghiselli F, Guerra D, Breton S, Passamonti M. 2013. A compar- Gamble T, Zarkower D. 2012. Sex determination. Curr Biol. ative analysis of mitochondrial ORFans: new clues on their origin and 22(8):R257–R262. role in species with doubly uniparental inheritance of mitochondria. Gangloff MM, Lenertz KK, Feminella JW. 2008. Parasitic mite and trematode Genome Biol Evol. 5(7):1408–1434. abundance are associated with reduced reproductive output and phys- Milani L, Ghiselli F, Maurizii MG, Passamonti M. 2011. Doubly uniparental iological condition of freshwater mussels. Hydrobiologia 610(1):25. inheritance of mitochondria as a model system for studying germ line Ghiselli F, et al. 2012. De novo assembly of the manila clam Ruditapes formation. PLoS One 6(11):e28194. philippinarum transcriptome provides new insights into expression Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 589 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Capt et al. GBE Milani L, Ghiselli F, Passamonti M. 2016. Mitochondrial selﬁsh elements Soneson C, Delorenzi M. 2013. A comparison of methods for differential and the evolution of biological novelties. Curr Zool. 62(6):687–697. expression analysis of RNA-seq data. BMC Bioinformatics 14:91. Milani L, Ghiselli F, Pecci A, Maurizii MG, Passamonti M. 2015. The ex- Sutovsky P, et al. 1999. Ubiquitin tag for sperm mitochondria. Nature pression of a novel mitochondrially-encoded gene in gonadic precur- 402(6760):371–372. sors may drive paternal inheritance of mitochondria. PLoS One Teaniniuraitemoana V, et al. 2014. Gonad transcriptome analysis of pearl 10(9):e0137468. oyster Pinctada margaritifera: identiﬁcation of potential sex differenti- Mitchell A, Guerra D, Stewart D, Breton S. 2016. In silico analyses of ation and sex determining genes. BMC Genomics 15:491. mitochondrial ORFans in freshwater mussels (Bivalvia: Unionoida) pro- The UniProt Consortium 2017. UniProt: the universal protein knowledge- vide a framework for future studies of their origin and function. BMC base. Nucleic Acids Res. 45(D1):D158–D169. Genomics 17:597. Tong Y, et al. 2015. Transcriptomics analysis of Crassostrea hongkongensis Moy GW, Springer SA, Adams SL, Swanson WJ, Vacquier VD. 2008. for the discovery of reproduction-related genes. PLoS One Extraordinary intraspeciﬁc diversity in oyster sperm bindin. Proc Natl 10(8):e0134280. Acad Sci U S A. 105:1993–1998. Uhlenhaut NH, Treier M. 2006. Foxl2 function in ovarian development. Mu ¨ ller T, et al. 2015. Factors affecting trematode infection rates in fresh- Mol Genet Metab. 88(3):225–234. water mussels. Hydrobiologia 742(1):59–70. Umen JG, Goodenough UW. 2001. Chloroplast DNA methylation and Negri I, et al. 2009. Unravelling the Wolbachia evolutionary role: the inheritance in Chlamydomonas. Genes Dev. 15(19):2585–2597. reprogramming of the host genomic imprinting. Proc R Soc B Biol Wallis MC, Waters PD, Graves JAM. 2008. Sex determination in Sci. 276(1666):2485–2491. mammals—before and after the evolution of SRY. Cell Mol Life Sci. Passamonti M, Ghiselli F. 2009. Doubly uniparental inheritance: two mi- 65(20):3182–3195. tochondrial genomes, one precious model for organelle DNA inheri- Wang X, et al. 2017. Transcriptome analysis of the freshwater pearl mussel tance and evolution. DNA Cell Biol. 28(2):79–89. (Cristaria plicata) mantle unravels genes involved in the formation of Patnaik BB, et al. 2016. Sequencing, de novo assembly, and annotation of shell and pearl. Mol Genet Genomics 292(2):343–352. the transcriptome of the endangered freshwater pearl bivalve, Cristaria Wijayawardena BK, Minchella DJ, DeWoody JA. 2016. The inﬂuence of plicata, provides novel insights into functional genes and marker dis- trematode parasite burden on gene expression in a mammalian host. covery. PLoS One 11(2):e0148622. BMC Genomics 17(1):600. Peng J, et al. 2015. Gonadal transcriptomic analysis and differentially Yarra T, Gharbi K, Blaxter M, Peck LS, Clark MS. 2016. Characterization of expressed genes in the testis and ovary of the Paciﬁc white shrimp the mantle transcriptome in bivalves: Pecten maximus, Mytilus edulis (Litopenaeus vannamei). BMC Genomics 16:1006. and Crassostrea gigas. Mar Genomics 27:9–15. Perlman SJ, Hodson CN, Hamilton PT, Opit GP, Gowen BE. 2015. Maternal Young MD, Wakeﬁeld MJ, Smyth GK, Oshlack A. 2010. Gene ontology transmission, sex ratio distortion, and mitochondria. Proc Natl Acad Sci analysis for RNA-seq: accounting for selection bias. Genome Biol. U S A. 112(33):10162–10168. 11(2):R14. Pozzi A, Plazzi F, Milani L, Ghiselli F, Passamonti M. 2017. SmithRNAs: Yu F-F, Wang M-F, Zhou L, Gui J-F, Yu X-Y. 2011. Molecular cloning and could mitochondria “bend” nuclear regulation? Mol Biol Evol. expression characterization of Dmrt2 in akoya pearl oysters, Pinctada 34(8):1960–1973. martensii. J Shellﬁsh Res. 30(2):247–254. Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor Zhang G, et al. 2012. The oyster genome reveals stress adaptation and package for differential expression analysis of digital gene expression complexity of shell formation. Nature 490(7418):49–54. data. Bioinformatics 26(1):139–140. Zhang ZH, et al. 2014. A comparative study of techniques for differential Robinson MD, Oshlack A. 2010. A scaling normalization method for dif- expression analysis of RNA-seq data. PLoS One 9:e1360 doi: 10.1371/ ferential expression analysis of RNA-seq data. Genome Biol. 11(3):R25. journal.pone.0103207 Sager R, Grabowy C. 1983. Differential methylation of chloroplast DNA Zhang N, Xu F, Guo X. 2014. Genomic analysis of the Paciﬁc oyster regulates maternal inheritance in a methylated mutant of (Crassostrea gigas) reveals possible conservation of vertebrate sex de- Chlamydomonas. Proc Natl Acad Sci U S A. 80(10):3025–3029. termination in a mollusc. G3 (Bethesda) 4:2207–2217. Santerre C, Sourdaine P, Adeline B, Martinez A-S. 2014. Cg-SoxE and Cg- Zhao X, Yu H, Kong L, Liu S, Li Q. 2014. Comparative transcriptome b-catenin, two new potential actors of the sex-determining pathway analysis of two oysters, Crassostrea gigas and Crassostrea hongkon- in a hermaphrodite lophotrochozoan, the Paciﬁc oyster Crassostrea gensis provides insights into adaptation to hypo-osmotic conditions. gigas. Comp Biochem Physiol Part A Mol Integr Physiol. 167:68–76. PLoS One 9(11):e111915. Sayadi A, Immonen E, Bayram H, Arnqvist G, Falabella P. 2016. The de Zieritz A, Aldridge DC. 2011. Sexual, habitat-constrained and parasite- novo transcriptome and its functional annotation in the seed beetle induced dimorphism in the shell of a freshwater mussel (Anodonta Callosobruchus maculatus. PLoS One 11(7):e0158565. anatina, Unionidae). J Morphol. 272(11):1365–1375. Shi J, Hong Y, Sheng J, Peng K, Wang J. 2015. De novo transcriptome Zouros E. 2013. biparental inheritance through uniparental transmission: sequencing to identify the sex-determination genes in Hyriopsis schle- the doubly uniparental inheritance (DUI) of mitochondrial DNA. Evol gelii. Biosci Biotechnol Biochem. 0:1–9. Biol. 40(1):1–31. Sim~ ao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31(19): 3210–3212. Associate editor: Bill Martin 590 Genome Biol. Evol. 10(2):577–590 doi:10.1093/gbe/evy019 Advance Access publication January 19, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/577/4817509 by Ed 'DeepDyve' Gillespie user on 16 March 2018
Genome Biology and Evolution – Oxford University Press
Published: Feb 1, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud