Pursuing the quest for better understanding the taxonomic distribution of the system of doubly uniparental inheritance of mtDNA

Specimen’s collection

Mature specimens of Yoldia hyperborea(Gould, 1841) were collected in the Baffin Sea (76°20^′50N, 77°35^′86W) in August 2013. Scrobicularia plana(Da Costa, 1778) samples were directly sent from the French National Museum of National History to our laboratory. All specimens were conserved in 95% ethanol. To identify sex-biased heteroplasmy, each individual was sexed by inspecting the gonads under a light microscope (100X) for the presence of eggs or sperm, and only individuals unambiguously sexed were kept for the present study. Dissections were carried on each individual to obtain somatic tissues (i.e., gills) and female or male gonad for DNA extractions (see below). The number of specimens analyzed for both species include 7 males and 7 females.

DNA extraction, polymerase chain reaction amplification and sequencing

Total genomic DNA was extracted separately from gonad tissue and from gill tissue with a Qiagen DNeasy Blood & Tissue Kit (QIAGEN Inc., Valencia, CA, USA) using the animal tissue protocol. The quality and quantity of DNA, respectively, were assessed by electrophoresis on 1% agarose gels and with a BioDrop µLITE spectrophotometer. Before PCR amplifications, all samples were treated using OneStep™ PCR Inhibitor Removal Kit (Zymo Research, Irvine, CA) according to the manufacturer’s protocol. For both species, partial sequence amplification of cytochrome oxidase subunit 1 (cox1) and large subunit ribosomal RNA (rrnL or 16S) were carried out in 50 µl volumes comprising 5.0 µl 10X Taq buffer, 1.0 µl dNTP mix (10mM), 2.0 µl of each forward and reverse primer (10 µM; LCO1490 and HCO2198 for cox1(Folmer et al., 1994), and 16Sar and 16Sbr for rrnL(Palumbi et al., 1991)), 0.25 µl Taq DNA Polymerase (5U/µl; Bio Basic Inc., Markham, Ontario), 2 µl of DNA extract (100 ng/ul), and ddH2O up to 50 µl. Reactions were performed on a TProfessional Basic Thermocycler with the following PCR amplification conditions: initial denaturation at 95°C for 2 min, followed by 35 cycles of 95°C for 20 sec, 44°C for 40 sec for cox1 and temperature gradient 40–60°C for 40 sec for rrnL, and 72°C for 40 sec, followed by a final extension step at 72°C for 5 min. The universal primers 16Sar and 16Sbr failed to amplify the rrnL sequence in both S. plana and Y. hyperborea. Resulting PCR products for cox1 were visualized on 1% agarose gels under UV light with SYBR green dye (Life Technologies), and purified with the Qiagen QIAquick PCR Purification Kit according to the manufacturer’s protocol. The purified PCR products were sequenced at the Genome Quebec Innovation Centre (McGill University), using the Applied Biosystem’s 3730xl DNA Analyzer technology.

DNA cloning and sequencing

Examination of chromatograms revealed the presence of multiple sequencing peaks only in male gonad tissues of S. plana and Y. hyperborea, suggesting co-amplification of two different mtDNA types. The amplified products of male gonads were thus cloned using the PGEM-T Easy vector (Promega, Madison, WI, USA) to confirm the presence of F and M genomes. Ten recombinant clones, for each species, were sent to the Genome Quebec Innovation Centre to be sequenced on both strands using the primers pUC20 (5^′-GTTTTCCCAGTCACGAC-3^′) and pUC2 (5^′-GAGCGGATAACAATTTCAC-3^′).

Sequence analysis

Cox1 sequences were edited and aligned using MEGA 6 (version6.06; Tamura et al., 2013). Amino acid sequences were deduced using the invertebrate mitochondrial genetic code. Calculations of nucleotide and amino acid p-distances were performed with MEGA 6 (with 1,000 bootstrap replicates) (version6.06; Tamura et al., 2013).

Following a similar approach than Dégletagne, Abele & Held (2016) to look for the presence of two intraspecific “F and M” clades in S. plana and Y. hyperborea, maximum likelihood (using RAxML version 8.2.8; Stamakis, 2014) with bootstrap analyses (1,000 replicates) and Bayesian phylogenies (using MrBayes v3.2.6; Huelsenbeck & Ronquist, 2001; Huelsenbeck & Ronquist, 2005; Ronquist et al., 2012) were performed on cox1 nucleotide sequences of both species with Soletellina virescens(Bivalvia, Veneridae, Genbank accession number: JN859944) and Yoldia eightsii(Bivalvia, Nuculanida, Genbank accession number: KJ571167) as outgroups for S. plana and Y. hyperborea, respectively (i.e., closest sequences according to BLAST search). Bayesian Information Criterion (BIC) (Schwarz, 1978) implemented in PartitionFinder (v1.1.1; Lanfear et al., 2012) was used to estimates the best-fitting models of evolution. Figtree (v1.4.2; Morariu et al., 2008) was used to edit the phylogenetic trees.

S. plana and Y hyperborea were included in an expanded phylogenetic analyses designed to verify molecular relationships among DUI species in general. Maximum likelihood (ML), maximum parsimony (MP) and Bayesian phylogenies were thus performed on F and M cox1 nucleotide and amino acid sequences from all DUI species known to date and Octopus vulgaris(Mollusca: Octopoda) and Aplysia californica(Mollusca: Gastropoda) were used as outgroup taxa. F and M cox1 sequences other than those obtained in the present study for S. plana and Y. hyperborea were retrieved from Genbank; the complete phylogenetic dataset is shown in Table 1. Cox1 sequences were aligned using MEGA 6 (version6.06; Tamura et al., 2013) and the best-fitting models of DNA evolution were selected using PartitionFinder (v1.1.1; Lanfear et al., 2012) according to BIC values (Schwartz,1978). Best models were applied whenever possible. Data were partitioned according to nucleotide position and gaps were treated as missing data.

ML analyses were conducted with RAxML (version 8.2.8; Stamatakis, 2014). A non-parametric bootstrap (Felsenstein, 1985) analysis was performed, using 1,000 bootstrap replicates and 20 ML searches, to assess nodal support for both trees. Outgroups were set to be paraphyletic to the monophyletic ingroup. MP analyses were carried out using PAUP software (v 4.0a147; Swofford, 2001). To optimize the chance of having the best topology, 100 random stepwise additions under tree-bisection reconnection branch swapping were applied (Bogan & Hoeh, 2000). Reliability of the internal nodes was evaluated by 1,000 pseudoreplicates using the heuristic search algorithm. Bayesian analyses were conducted using MrBayes (v3.2.6; Huelsenbeck & Ronquist, 2001; Huelsenbeck & Ronquist, 2005; Ronquist et al., 2012). Each analysis consisted of two independent runs of 4 MC³ chains that were run for 10,000,000 generations. Convergence was estimated through the log likelihood value of trees, potential scale reduction factor (PSRF) and standard deviation of average split frequencies sampled every 1,000 generations (Gelman & Rubin, 1992). Trees were sampled every 100 generations and a majority-rule consensus tree was computed after discarding the first 25% as burn-in. Fidelity of the topology was evaluated with the posterior probabilities from the consensus tree. All phylogenetic trees were edited for easier readability using FigTree (v1.4.2; Morariu et al., 2008).

All the alignments used for phylogenetic reconstruction are available here (https://dx.doi.org/10.6084/m9.figshare.3798789.v1).

Genetic distances

For this study, two new species were tested for the presence of DUI and 42 sequences were examined in total: the number of sequences and haplotypes for each species are listed in Table 2. Different haplotypes (i.e., for the same tissue among females or males, and thus non related to DUI) were found for both species (Table 2). All mtDNA sequences are available via GenBank under accession numbers KX447420, KX447421, KX447422, KX447423, KX447424, KX447425, KX447426, KX447427, KX447428. Sequences with the same haplotype were deposited only once.

Intragroup (female sequences and male sequences, respectively) and intergroup (female versus male sequences) nucleotide and amino acid p-distances are shown for both species in Table 3. The p-distances between female and male sequences within each species are significantly larger than the within group p-distances. Specifically, for Y. hyperborea cox1 sequences, intragroup p-distances are low, i.e., 0.0014 for female sequences and zero for male sequences (i.e., male sequences are identical), with standard error of ±0.0008 for female sequences, whereas the between group p-distance is considerably higher with a value of 0.0596 ± 0.0079 (amino acid p-distance is 0.0454 ± 0.0124). The same observation can be made for S. plana cox1 sequences: within group p-distances are 0.0067 ± 0.0034 and 0.0020 ± 0.0014, for female and male sequences, respectively, whereas the between group p-distance value is 0.0965  ± 0.0074 (amino acid p-distance is 0.0659 ± 0.0100).

Phylogenetic analyses for Scrobicularia plana and Yoldia hyperborea

Phylogenetic analyses of partially sequenced S. plana cox1 and Y. hyperborea cox1 sequences were conducted using the HKY + G model (Hasegawa, Kishino & Yano, 1985). Based on high bootstraps and posterior probability values, female and male haplotypes were clustered into two well-supported clades for each species (Figs. S1 and S2).

Taxonomic distribution of DUI

Genetic distances and phylogenetic analyses provide evidence for sex-associated mitochondrial heteroplasmy in Scrobicularia plana(Bivalvia: Semelidae) and Yoldia hyperborea(Bivalvia: Yoldiidae), raising the total of bivalve families in which DUI has been discovered to 12 (Hoeh, Stewart & Guttman, 2002; Theologidis et al., 2008; Boyle & Etter, 2013; Plazzi, 2015; Dégletagne, Abele & Held, 2016).

The list of the 103 bivalve species in which DUI has been found to date is presented in Table 4. Evidence for DUI is reported in the literature for 92 species (and for two species in the present study). For the other nine species, the evidence is based on sequences derived from male and female gonads and retrieved from GenBank (see Tables 1 and 4). These sequences all show strong nucleotide divergence between mitochondrial gene sequences obtained from male and female individuals, with F to M type p-distances ranging from 10% to 30%.

Phylogenetic analyses

A total of 608 cox1 nucleotide sequences were aligned for a total of 114 operational taxonomic units (OTUs) for phylogenetic analyses. GTR + I + G (Tavaré, 1986) was selected as the best-fitting model of evolution for nucleotides (Table S1). The majority rule nucleotide-based BI tree was favored as our best tree obtained for this study (Figs. 1–4). It shows strong topology similarity with both MP and ML trees (Figs. S3 and S4).

Main features of the BI tree (Figs. 1–4) are as follows: (1) pteriomorph (Mytiloida) + Yoldia hyperborea (PP = 0,80) and heterodont bivalves (Veneroida) (PP = 0,99) are reciprocally monophyletic with palaeoheterodont (Unionoida) bivalves being the sister group to these clades (PP = 0,86); (2) the order Veneroida is well resolved with 3 nodes (PP = 1.00) separating each superfamily: Veneroidea, Tellinoidea and Mactroidea; (3) in both F and M clusters of the order Unionoida, the superfamily Hyrioidea represented by H. menziesi is a positioned as a well separated sister group (PP = 1.00) to the Unionoidea superfamily (Note: within the Unionoidea, the separation between the two families Margaritiferidae and Unionidae is also apparent and well supported only in the M cluster (PP = 0.93)); (4) the F and M clades are reciprocally monophyletic only in unionoids and Mytilus spp., i.e., the F sequences of different species cluster together as do the M sequences, all the other species exhibit a phylogenetic pattern where F and M sequences clusters are distinct from one another but are nonetheless sister groups; and (5) branch lengths indicate a higher substitution rate for the M genomes relative to that of the F genomes for almost all species.

Evidence for DUI in Scrobicularia plana and Yoldia hyperborea

Both p-distance and phylogenetic analyses indicate the coexistence of sex-linked F and M mitochondrial lineages in S. plana and Y. hyperborea(Figs. 1–4; Table 3 and Figs. S1 and S2). As mentioned earlier, the strategy of searching for sex-biased heteroplasmy of mitochondrial gene sequences as a means for detecting DUI has been successfully employed in other bivalve species (e.g., Boyle & Etter, 2013; Vargas et al., 2015; Dégletagne, Abele & Held, 2016). DUI can also be detected by in vivo localization of male mitochondria in embryos: an aggregate pattern of M-type mitochondria only into the male germline during early embryonic stages is typical of DUI species (Cao, Kenchington & Zouros, 2004; Obata & Komaru, 2005; Cogswell, Kenchington & Zouros, 2006; Milani, Ghiselli & Passamonti, 2012). The percentage of nucleotide divergence between the F and M cox1 sequences for both of these species, i.e., 33.4% for S. plana and 13% for Y. hyperborea, is within the range of what has been found for other bivalves with DUI (e.g., 8% in the veneroid Artica islandica(Dégletagne, Abele & Held, 2016), 17% in the nuculanoid Ledella sublevis(Boyle & Etter, 2013); 24% in the mytiloid Mytilus edulis(Breton et al., 2006),and 50% in the unionoid Inversidens japanensis(Doucet-Beaupré et al., 2010)).

The discovery of DUI in Y. hyperborea provides the first example of this unusual system of mitochondrial DNA transmission in the family Yoldiidae, and the third case for the order Nuculanoida (subclass Protobranchia; Boyle & Etter, 2013). Protobranchia being the most basal lineage within the Bivalvia (Giribet & Wheeler, 2002; Smith et al., 2011), this result suggests that the origin of DUI is older than the birth of the autolamellibranchiata (Theologidis et al., 2008; Doucet-Beaupré et al., 2010; Boyle & Etter, 2013). However, of the three protobranch orders, only the Nuculanoida, has been reported to exhibit DUI (Boyle & Etter, 2013), and the phylogenetic status of this order is still being questioned. Indeed, recent phylogenetic studies suggest that the Nuculanoida is not a member of the basal protobranch group, which includes Nuculoida and Solemyoida, but instead is associated with the Pteriomorpha (Wilson, Rouse & Giribet, 2010 (based on 3 nuclear and 2 mitochondrial genes); Plazzi et al., 2011) (based on 4 mitochondrial genes); Breton et al., unpublished (based on complete mitochondrial genomes and 3 nuclear genes combined with an extensive morphological dataset)), whereas other recent studies based on four nuclear genes (Sharma et al., 2012) and phylogenomic data (Smith et al., 2011; González et al., 2015) supported the monophyly of Nuculanoida + Opponobranchia (Nuculoida and Solemyoida; Giribet, 2008). The presence of DUI in protobranchs thus remains an open question until the publication of a well-supported and robust phylogeny of bivalves showing the monophyly of the traditional clade Protobranchia (i.e., Solemyoida + Nuculoida + Nuculanoida) and/or until the discovery of DUI in nuculoid or solemyoid bivalves.

The peppery furrow shell Scrobicularia plana belongs to the order Veneroida. It is the first reported species with DUI from the family Semelidae, raising the total number of veneroid families in which DUI has been discovered to six (Theologidis et al., 2008; Plazzi, 2015; Dégletagne, Abele & Held, 2016; present study). Apart from providing new insights into the taxonomic distribution of DUI, S. plana may play a key role for better understanding the hypothesized role of DUI in sex determination (e.g., Breton et al., 2007; Breton et al., 2011a; Breton et al., 2014; Breton & Stewart, 2015; Mitchell et al., 2016). Indeed, an “intersex” condition, i.e., the appearance of oocytes in male gonads following endocrine disruption, has been reported in S. plana and is associated with differentially expressed mitochondrial transcripts in males exhibiting intersex compared to “normal” males (Chesman & Langston, 2006). Specifically, using a suppressive subtractive hybridization approach, Ciocan et al. (2012) were able to determine that several mitochondrial mRNA transcripts were down-regulated in clam intersex samples (i.e., cox1, cytb, nad1, nad2, nad3, nad4). Interestingly, we observed that the down-regulated cox1 sequence identified by Ciocan et al. (2012) was identical to the male cox1 sequences from our study, indicating that the down-regulation of male mitochondrial sequences is associated with the appearance of female characteristics in male gonads in this species. These results provide more evidence for a link between DUI and sex determination or differentiation. It is noteworthy that the intersex has been shown to be a widespread phenomenon in bivalves, including in species with DUI (e.g., R. philippinarum(Lee et al., 2010) and M. galloprovincialis (Ortiz-Zarragoitia & Cajaraville, 2010)).

Taxonomic distribution of DUI in bivalves: an update

Including the two species in the present study, DUI has been reported to date in 104 bivalve species belonging to four subclasses (Heterodonta, Palaeoheterodonta, Pteriomorphia, Protobranchia (but see above comments regarding the questionable inclusion of Nuculanoida within the Protobranchia)), four orders (Mytiloida, Nuculanoida, Unionoida, Veneroida), nine superfamilies (Arcticoidea, Hyrioidea, Mactroidea, Mytiloidea, Nuculanoidea, Solenoidea, Tellinoidea, Unionoidea, Veneroidea) and twelve families (Arcticidae, Donacidae, Hyriidae, Mactridae Margaritiferidae, Mytilidae, Nuculanidae, Semelidae, Solenidae, Unionidae Veneridae, Yoldiidae) (Table 4 and Figs. 1–4). However, DUI is certainly more widespread in the Bivalvia given that its detection remains difficult; the higher rate of molecular evolution of M type mitochondrial genomes may make it less likely that “universal” mitochondrial primers will anneal to and amplify the M type (Theologidis et al., 2008; Zouros, 2013). In addition, the process of mitochondrial genome “masculinization,” i.e., when an F genome invades the male route of transmission, can also make the paternally-transmitted genome almost indistinguishable from the maternally-transmitted one (Stewart et al., 2009; Theologidis et al., 2008; Zouros, 2013). Additional studies of bivalves and other mollusc species will significantly contribute to better understanding the taxonomic distribution of the system of doubly uniparental inheritance of mtDNA.

Phylogenetic analyses and the origin of DUI

As for other veneroid and nuculanoid species (e.g., Theologidis et al., 2008), the sex-linked mtDNA sequences of S. plana and Y. hyperborean, respectively, exhibit a phylogenetic pattern in which the F and M mtDNA sequences are different from one another but yet cluster together in a monophyletic group (Figs. 1–4). Such a pattern can also be seen in the order Mytiloida (Figs. 1–4), except for the Mytilus species complex, which is in agreement with previous studies (e.g., Doucet-Beaupré et al., 2010). On the other hand, the observed F/M phylogeny of unionoids contrasts with the patterns observed in the mytiloids, nuculanoids and veneroids. In unionoids, all of the F sequences cluster together and all of the M genomes cluster together such that the F sequences form a monophyletic clade and the M sequences form a monophyletic clade. Similar results have previously been obtained and have suggested that the M and F lineages in the order Unionoida have been distinct and maintained for at least 200 million years (e.g., Doucet-Beaupré et al., 2010). The observed phylogenetic pattern of the unionoids could be related to the absence of a masculinization event in this group for over 200 million years (Hoeh, Stewart & Guttman, 2002). One hypothesis explaining why F-to-M masculinization events do not occur in freshwater mussels involves the cox2 extension present only in the M genome (Curole & Kocher, 2002). If this extension is essential for the function or transmission of the M genome (Curole & Kocher, 2002; Chakrabarti et al., 2006; Chakrabarti et al., 2007) then masculinization events are no longer possible in this group, unless F genomes are able to obtain, by recombination, those male-specific elements necessary for being paternally transmitted (Stewart et al., 2009; Doucet-Beaupré et al., 2010). The M genome is considered as a “nearly selfish” element in the sense that it does fulfill its function only if this work cannot be done by the F genome. In other words, paternally-transmitted mt genomes only perform male specific functions. This restricted functionality of the M mtDNA to male gonad tissue is one factor that may explain why the M genome usually evolves faster than the F genome in species with DUI, i.e., because of relaxed selective constraints (Stewart et al., 1996; Zouros, 2000; Passamonti, Boore & Scali, 2003). The other hypotheses that have been proposed to explain the higher rate of evolution of the M genome include (i) a higher number of male germ line mitotic divisions preceding gametogenesis compared to the female germ line (Selwood, 1968; Zwaan & Mathieu, 1992; Stewart et al., 1996), (ii) a greater degree of oxidative damage in metabolically active sperm relative to eggs (Stewart et al., 1996; Zouros, 2013) and/or (iii) a smaller effective population size of male mitochondria compared to female mitochondria (Stewart et al., 1996; Zouros, 2013).

As for other recent morphological and molecular studies based on mtDNA and nuclear sequences (Giribet & Wheeler, 2002; Giribet & Distel, 2003; Bieler & Mikkelsen, 2006; Plazzi & Passamonti, 2010; Plazzi et al., 2011; S Breton et al., 2016, unpublished data), our results support a more derived Nuculanoida clustering with Pteriomorphia instead of a basal position of Nuculanoida + Opponobranchia (=Protobranchia). Therefore, the question about the origin of DUI in the branch leading to the Autolamellibranchia about 460 Mya (Little & Vrijenhoek, 2003) or much earlier, perhaps in the early Cambrian (Boyle & Etter, 2013), still remains open. The taxonomic position of Y. hyperborea and other members of the order Nuculanoida should be tested using an expanded data set. A robust bivalve phylogeny, as well as a much more accurate understanding of the taxonomic distribution of DUI, is needed to allow a rigorous evaluation of a single vs. multiple origins of DUI.

Figures 1–4 demonstrate one split between M and F types at the base of the Unionoida, and one split for each species exhibiting a taxon-specific pattern. Taken at face value, this pattern implies a large number of independent origins of DUI. Given the rare molecular and developmental complexity associated to the DUI system, this hypothesis seems unlikely but cannot be completely rejected yet (Zouros, 2013). The opposite hypothesis of a single origin of DUI can only be true if associated with masculinization events along each branch of the phylogenetic tree of the Bivalvia where F and M types are each other’s closest relatives. Such masculinization events have already been clearly demonstrated in Mytilus (reviewed in Stewart et al., 2009; Zouros, 2013), and evidence suggests that mitochondrial recombination and acquisition of key elements of the evolutionarily older M mt genome (i.e., sperm transmission elements) are necessary for a F genome to be transmitted via sperm (e.g., Stewart et al., 2009; Zouros, 2013; Kyriakou et al., 2015). Because masculinization events restore nucleotide divergence between F and M mtDNAs to zero for most of the genes in each of the sex-associated genome (except for key sperm transmission elements), this phenomenon could explain the F/M phylogenetic patterns of mytiloids, nuculanoids and veneroids. As a consequence of the many similarities found among the distantly-related DUI species (e.g., sex ratio bias, mitochondria’s behavior in the newly formed zygotes, rates of evolution of the two genomes), and because of the complexity of the DUI system, one could thus favor the hypothesis of a single origin of DUI with repeated masculinization events. However, recent in silico analyses of the novel mtDNA-encoded protein-coding genes, i.e., ORFan genes with no known homologous, which have been found in species with DUI support either a viral or a mitochondrial origin for these genes that are most probably involved in the DUI mechanism (Milani et al., 2013; Milani et al., 2014; Mitchell et al., 2016; Milani, Ghiselli & Passamonti, in press). These results suggest the possibility of DUI systems with elements of different sources/origins and different mechanisms of action in the distantly-related DUI taxa (i.e., DUI could be achieved by different modifications of strictly maternal inheritance of mitochondria), and this scenario would best fit multiple origins of DUI and the necessary factors.

To conclude, our study presents evidence for the existence of DUI in the nuculanoid species Yoldia hyperborea and the veneroid species Scrobicularia plana. Because the taxonomic position of Y. hyperborea and its order Nuculanoida has been debated over the years and still remain uncertain, the question about the origin of DUI during the Cambrian or before is still unresolved. A much more accurate understanding of the taxonomic distribution of DUI across the Bivalvia appears to be a priority to help confirming the single or multiple origins of this unusual system of mitochondrial heredity. Testing for its presence in all bivalve superfamilies, especially the one not studied yet (e.g., Solemyoidea, Lucinoidea, Carditoidea), is the only way to make a definitive statement.