Mitochondrial Architecture Rearrangements Produce Asymmetrical Nonadaptive Mutational Pressures That Subvert the Phylogenetic Reconstruction in Isopoda.

The phylogeny of Isopoda, a speciose order of crustaceans, remains unresolved, with different data sets (morphological, nuclear, mitochondrial) often producing starkly incongruent phylogenetic hypotheses. We hypothesized that extreme diversity in their life histories might be causing compositional heterogeneity/heterotachy in their mitochondrial genomes, and compromising the phylogenetic reconstruction. We tested the effects of different data sets (mitochondrial, nuclear, nucleotides, amino acids, concatenated genes, individual genes, gene orders), phylogenetic algorithms (assuming data homogeneity, heterogeneity, and heterotachy), and partitioning; and found that almost all of them produced unique topologies. As we also found that mitogenomes of Asellota and two Cymothoida families (Cymothoidae and Corallanidae) possess inversed base (GC) skew patterns in comparison to other isopods, we concluded that inverted skews cause long-branch attraction phylogenetic artifacts between these taxa. These asymmetrical skews are most likely driven by multiple independent inversions of origin of replication (i.e., nonadaptive mutational pressures). Although the PhyloBayes CAT-GTR algorithm managed to attenuate some of these artifacts (and outperform partitioning), mitochondrial data have limited applicability for reconstructing the phylogeny of Isopoda. Regardless of this, our analyses allowed us to propose solutions to some unresolved phylogenetic debates, and support Asellota are the most likely candidate for the basal isopod branch. As our findings show that architectural rearrangements might produce major compositional biases even on relatively short evolutionary timescales, the implications are that proving the suitability of data via composition skew analyses should be a prerequisite for every study that aims to use mitochondrial data for phylogenetic reconstruction, even among closely related taxa.

Introduction

Significant taxonomic and phylogenetic uncertainty permeates the entire order of Isopoda, a highly speciose (>10,000) order of crustaceans (class Malacostraca) that exhibit a remarkable diversity in their life histories and occupy almost all habitats on the planet Earth. The traditional morphology-based taxonomic classification and identification of isopods is further (aside from the speciosity) hampered by great intraspecific morphological variation, sexual dimorphism, sequential hermaphroditism, relatively flexible host preference, and global distribution of many species (Wetzer 2002; Wilson 2008; Lins et al. 2012, 2017; Joca et al. 2015; Shen et al. 2017; Rudy et al. 2018). However, molecular data also appear to be an unreliable tool for the task, as different data sets (mitochondrial genes, mitochondrial genomes, nuclear genes, combined mitonuclear data) often produce very different topologies (Brusca 1981; Wetzer 2002; Wilson 2009; Kilpert et al. 2012; Poore and Bruce 2012; Martin et al. 2016; Hata et al. 2017; Lins et al. 2017; Yu et al. 2018; Zou et al. 2018). As a result, even the identity of the basal isopod clade (defined as the sister-clade to all other isopod lineages; Krell and Cranston 2004) remains debated. Traditionally (morphology and single gene-based studies), Phreatoicidea was regarded as the basal clade (Wetzer 2002; Wilson 2009; Kilpert et al. 2012). However, some studies resolved Phreatoicidea+Aselotta at the base (Wilson 1999, 2009; Shen et al. 2017; Yu et al. 2018), one study found Limnoriidea (Lins et al. 2017) at the base, whereas a few studies even resolved parasitic Cymothoidae and Corallanidae (suborder Cymothoida) at the base (Wilson 2009; Lins et al. 2012, 2017; Hua et al. 2018; Zou et al. 2018). As the suborder Cymothoida was traditionally regarded as the most derived isopod clade (Brusca 1981; Wetzer 2002; Wilson 2009; Kilpert et al. 2012), these alternative hypotheses cannot be described as a minor topological instability. The monophyly of the Cymothoida is generally supported by the morphological data, but rejected by the molecular data (Brandt and Poore 2003; Wilson 2009; Kilpert et al. 2012; Lins et al. 2017; Shen et al. 2017; Hua et al. 2018; Yu et al. 2018; Zou et al. 2018). Among a number of other unresolved phylogenetic issues permeating this order are the monophyly of the suborder Oniscidea (supported by morphology, sometimes rejected by the molecular data) and the existence of several “rogue” species/taxa, such as Ligia oceanica (Ligiidae), Eurydice pulchra (Cirolanidae), and Limnoria quadripunctata (Limnoriidae), whose positions in the isopod clade often vary among studies (Schmidt 2008; Wilson 2009; Kilpert et al. 2012; Wetzer et al. 2013; Lins et al. 2017; Shen et al. 2017; Yu et al. 2018).

Historically, variation found within gene sequences was typically considered to accumulate under a neutral equilibrium model, so commonly used phylogenetic reconstruction algorithms assume homogeneity in mutational rates. As this paradigm began to change during the last few decades (Wolff et al. 2014), this was accompanied by a growing amount of evidence that compositional heterogeneity can compromise phylogenetic reconstruction in some taxa. Therefore, the evolutionary models operating under that prerequisite may not be suitable for all phylogenetic studies (Kolaczkowski and Thornton 2004; Hassanin 2006; Phillips et al. 2006; Lartillot et al. 2007; Sheffield et al. 2009; Zhong et al. 2011; Morgan et al. 2013; Cameron 2014). However, the feud about the most suitable methodological approach to account for this heterogeneity remains unresolved, with most prominent contenders currently being the CAT models (Feuda et al. 2017), and partitioning schemes, that is, different evolutionary models assigned to different character blocks, assuming homogeneity within each block (CAT models to not allow partitioned data) (Whelan and Halanych 2017). Although these two approaches account for rate heterogeneity across sites, they still assume that substitution rates for sites are constant across all included lineages. From the evolutionary perspective, this is not a likely scenario, since substitution rates are likely to be both site- and lineage-specific (Crotty et al. 2017). Indeed, heterotachy, variations in lineage-specific evolutionary rates over time (Lopez et al. 2002), is widespread in eukaryotes (Baele et al. 2006).

Mitochondrial genomes (mitogenomes) generally provide much higher phylogenetic resolution than traditionally used morphological and single gene-based molecular markers (Nie et al. 2018), so mitochondrial phylogenomics is increasingly used to tackle phylogenetic controversies (Cameron 2014; Der Sarkissian et al. 2015; Lan et al. 2017; Li et al. 2017; Bourguignon et al. 2018). Although the resolution of this approach is still limited by a very small number of available mitogenomes in isopods, overview of published studies shows that they also failed to produce results congruent with other approaches and to resolve the rogue taxa issues (Kilpert et al. 2012; Lins et al. 2017; Shen et al. 2017; Hua et al. 2018; Yu et al. 2018; Zou et al. 2018). The evolutionary history of Isopoda abounds in independent (major) life history innovations, such as free-living to parasitic lifestyle (Jones et al. 2008; Ketmaier et al. 2008; Poore and Bruce 2012; Hata et al. 2017), and radical habitat expansions (Lins et al. 2012), such as sea to freshwater, or even water to land (Wilson 2008; Broly et al. 2013; Hata et al. 2017). It has been proposed that the Cymothoida may have originated in deep seas, subsequently expanded to shallow seas, and then to brackish and freshwater (likely on several independent occasions) (Hata et al. 2017). Mitochondrial genes are central to the energy production via the oxidative phosphorylation (Gawryluk et al. 2016), and signals of adaptation to high altitude (Mishmar et al. 2003; Hassanin et al. 2009; Scott et al. 2011), deep-sea environment (Almeida et al. 2015), and shifts in physiological demands (Hassanin 2006; Botero-Castro et al. 2018), have been identified in mitogenomes of a range of animals. It is therefore highly likely that radical adaptations to life in different environments, from the anoxic environment of deep sea-inhabiting isopod species (Lins et al. 2012) to terrestrial species, would produce strikingly different evolutionary pressures on genomes of species, and result in disparate evolutionary rates of mitochondrial genes. In agreement with this hypothesis are uneven evolutionary rates (dN/dS) observed among isopod mitogenomes (Shen et al. 2017) and different mutational rates of protein-coding genes (PCGs) encoded on the majority (or plus) strand among different lineages of isopods (Lloyd et al. 2015). Conflicting phylogenetic signals among different mitochondrial regions have been reported in a number of metazoan groups, which indicates that different mitochondrial regions can accumulate substitutions in ways that are difficult to model, and thus produce biased estimates of phylogeny (Meiklejohn et al. 2014). There is evidence that this compositional heterogeneity may be comparatively highly pronounced in mitogenomes of some arthropod taxa (Hassanin 2006; Cameron 2014; Liu et al. 2017). Despite these limitations of molecular data (Ballard and Whitlock 2004; Hassanin et al. 2005; Rubinoff et al. 2005; Talavera et al. 2011; Grechko 2013; Edwards et al. 2016; Willis 2017), most previous phylogenetic studies of Isopoda ignored those limitations, or attempted to ameliorate them by using such strategies as combined data sets (mtDNA, nuclear DNA, morphology) (Wilson 2009; Lins et al. 2017), amino acid sequences (Kilpert et al. 2012; Lins et al. 2017), or applying different models to each codon position (Hata et al. 2017). However, none of those studies attempted to use algorithms designed specifically to account for compositional heterogeneity/tachy, nor studied this problem directly. We hypothesized that the extreme life history diversity of isopods might be causing pronounced compositional heterogeneity/heterotachy in their mitogenomes, and interfere with the reconstruction of their phylogeny. To test this hypothesis, we used a number of different data sets: mitochondrial DNA (single genes, genomes, nucleotides, amino acids, gene orders) and nuclear DNA (18S); and methodological approaches: data set partitioning, maximum likelihood (ML), Bayesian inference (BI), parsimony, PhyloBayes (PB) CAT-GTR model (heterogeneous), and GHOST (heterotachous).

Materials and Methods

Data Sets

As a majority of available isopod mitogenomes are incomplete, we were faced with a trade-off between the amount of data used in the analysis and the number of species used: after removing 6 (too few genes and duplicates) of the 27 available isopod mitogenomes (Oct, 2018), the data set comprised 8 complete and 13 partial sequences (table 1 and supplementary file S1, Supplementary Material online). Out of 11 valid isopod suborders, our data set covers 7. To attempt to resolve the debated issue of the basal Isopod clade with maximum resolution, we used a relatively large number of outgroups for phylogenetic analyses: a basal arthropod, Limulus polyphemus (Lavrov et al. 2000), and a number of nonisopod Malacostraca: three Decapoda, two Stomatopoda, two Amphipoda, one Mysida, and one Euphausiacea species (Kilpert and Podsiadlowski 2006; Wilson 2009). We conducted phylogenetic analyses using the following data sets: NUC—nucleotides of concatenated 13 PCGs and 2 rRNA genes (rrnL and rrnS), AAs—concatenated amino acid sequences of 13 PCGs, 15 single-gene data sets (13 PCGs + 2 rRNAs), and gene orders. We also tested the performance of data partitioning, by conducting the same analyses on both nonpartitioned and partitioned data sets.

Table 1

Taxonomy, Length (bp), Base Composition (%), and Skews of Mitogenomes Used in the Analysis

Species	Suborder	Family	Acc. No.	Length	A + T	A	C	AT Skew	GC Skew
Isopoda
Asellus aquaticus	Asellota	Asellidae	GU130252	13,639	61.9	31	21.3	0.002	−0.122
Cymothoa indica	Cymothoida	Cymothoidae	MH396438	14,475	63.8	36	26.1	0.129	−0.442
Ichthyoxenos japonensis	Cymothoida	Cymothoidae	MF419233	15,440	72.7	37	18.7	0.026	−0.375
Tachaea chinensis	Cymothoida	Corallanidae	MF419232	14,616	72.8	38	18.5	0.055	−0.354
Eurydice pulchra	Cymothoida	Cirolanidae	GU130253	13,055	55.9	27	17.7	−0.052	0.198
Bathynomus sp.	Cymothoida	Cirolanidae	KU057374	14,965	58.7	27	17	−0.093	0.175
Gyge ovalis	Cymothoida	Bopyridae	NC_037467	14,268	59.6	27	17.8	−0.093	0.118
Eophreatoicus sp.	Phreatoicidea	Amphisopidae	NC_013976	14,994	69.6	31	11.4	−0.104	0.25
Sphaeroma serratum	Sphaeromatidea	Sphaeromatidae	GU130256	13,467	54.4	25	17.8	−0.069	0.219
Limnoria quadripunctata	Limnoriidea	Limnoriidae	NC_024054	16,515	66.3	30	13.6	−0.104	0.188
Idotea balthica	Valvifera	Idoteidae	DQ442915	14,247	61	28	16.3	−0.076	0.163
Glyptonotus antarcticus	Valvifera	Chaetiliidae	GU130254	13,809	65.4	32	16.6	−0.033	0.038
Trachelipus rathkii	Oniscidea	Trachelipodidae	MF187612	14,080	67.3	33	13.3	−0.029	0.184
Porcellio dilatatus	Oniscidea	Porcellionidae	KX289582	14,103	65.6	31	13.4	−0.067	0.224
Cylisticus convexus	Oniscidea	Cylisticidae	KR013002	14,154	67.8	33	12.9	−0.035	0.194
Porcellionides pruinosus	Oniscidea	Porcellionidae	KX289584	14,078	60.5	28	14.9	−0.09	0.248
Armadillidium album	Oniscidea	Armadillidiidae	KX289585	13,812	69.7	33	12.5	−0.045	0.172
Armadillidium nasatum	Oniscidea	Armadillidiidae	MF187611	13,943	68.1	33	13.4	−0.043	0.16
Armadillidium vulgare	Oniscidea	Armadillidiidae	MF187614	13,932	71.5	34	11.8	−0.043	0.174
Armadillidium vulgare	Oniscidea	Armadillidiidae	MF187613	13,955	71.3	34	11.8	−0.039	0.179
Ligia oceanica	Oniscidea	Ligiidae	NC_008412	15,289	60.9	29	17	−0.041	0.134
Species	Order	Suborder	Acc. No.	Length	A+T	A	C	AT Skew	GC Skew
Nonisopod Malacostraca and Limulus polyphemus
Atergatis floridus	Decapoda	Pleocyemata	NC_037201	16,180	69.3	33	20.3	−0.036	−0.319
Penaeus vannamei	Decapoda	Dendrobranchiata	NC_009626	15,990	67.7	33	19.2	−0.026	−0.192
Typhlatya miravetensis	Decapoda	Pleocyemata	NC_036335	15,865	66.2	36	22.5	0.076	−0.332
Squilla mantis	Stomatopoda	Unipeltata	NC_006081	15,994	70.2	35	16.8	−0.001	−0.13
Lysiosquillina maculata	Stomatopoda	Unipeltata	NC_007443	16,325	63.9	33	21.4	0.026	−0.185
Metacrangonyx repens	Amphipoda	Senticaudata	NC_019653	14,355	76.9	38	11.7	−0.025	−0.014
Neomysis japonica	Mysida	NA	NC_027510	17,652	74.5	37	13.8	−0.021	−0.085
Eulimnogammarus cyaneus	Amphipoda	Senticaudata	NC_033360	14,370	67.6	33	20.3	−0.019	−0.251
Euphausia pacifica	Euphausiacea	NA	NC_016184	16,898	72	36	16	0.004	−0.145
Limulus polyphemus	Xiphosurida	NA	NC_003057	14,985	67.6	38	22.7	0.111	−0.399

To test a signal from the nuclear data, we used partial 18S rRNA gene (1,223 bp aligned data set; supplementary file S2, Supplementary Material online). As some of the species from the mitogenomic data set were not available, we made sure to include representatives of all isopod suborders in the mitogenomic data set (55 species in total; supplementary file S1, Supplementary Material online). To obtain a more comparable data set, we sequenced the 18S gene (≈1,930 bp) of three parasitic Cymothoida species from the mtDNA data set: Cymothoa indica (Cymothoidae; GenBank accession number MK079664), Ichthyoxenos japonensis (Cymothoidae; MK542857), and Tachaea chinensis (Corallanidae; MK542858). To improve the resolution for the “problematic” Cymothoidae, we sequenced an additional species: Asotana magnifica (MK542856).

Data Manipulation and Analyses

PhyloSuite (Zhang et al. 2018) was used to batch-download all selected molecular data from the GenBank, extract genomic features, translate genes into amino acid sequences, semiautomatically reannotate ambiguously annotated tRNA genes with the help of the ARWEN (Laslett and Canbäck 2008) output, automatically replace the GenBank taxonomy with the WoRMS database taxonomy, as the latter tends to be more up to date (Costello et al. 2013), generate comparative genome statistics tables, and conduct phylogenetic analyses (Flowchart mode) using a number of incorporated plug-in programs.

We used MAFFT (Katoh and Toh 2008; Katoh and Standley 2013) to align sequences: nucleotide and amino acid sequences of PCGs were aligned in batches (using codon and normal-alignment modes, respectively) with “–auto” strategy, whereas rRNA genes (including the 18S) were aligned using Q-INS-i algorithm, which takes secondary structure information into account. Gblocks (Castresana 2000; Talavera and Castresana 2007) was used to remove ambiguously aligned regions from the concatenated alignments (default PhyloSuite parameter settings). After concatenating the alignments with PhyloSuite (all alignments in supplementary file S2, Supplementary Material online), PartitionFinder2 (Lanfear et al. 2012) was used to find the best data partitioning scheme and to select the best-fit evolutionary models for each partition (Akaike Information Criterion; supplementary file S2, Supplementary Material online), whereas ModelFinder (Kalyaanamoorthy et al. 2017) was used to select models for nonpartitioned data sets (Bayesian information criterion). χ2 test for the homogeneity of character composition of aligned sequences was performed using IQ-TREE 1.6.8 (Trifinopoulos et al. 2016).

We tested the performance of two standard homogeneous models, ML, BI, both on nonpartitioned and partitioned data, and two nonstandard heterogeneous (CAT-GTR) and heterotachous (GHOST) models (these two require the input data to be nonpartitioned). We also tested the performance of Parsimony model implemented in PAUP* 4.0 on nonpartitioned data sets (heuristic searching, TBR branch swapping, and 500 random addition sequence replicates, bootstrap: 1,000 pseudoreplicate data sets) (Swofford 2002). BI analyses were conducted using MrBayes 3.2.6 (Ronquist et al. 2012) (as PhyloSuite plug-in), with default settings until the stationarity was reached (3.8–8.3 × 106 generations depending on the data set, stationarity = average SD of split frequencies <0.01, estimated sample size >200, PSRF index ≈1). ML analyses were mostly carried out using a ML+rapid bootstrap algorithm with 1,000 replicates in RAxML (Stamatakis 2014), with the exception of ML analyses of 15 single-gene data sets and 18S, which were conducted in PhyloSuite batch mode using IQ-TREE plug-in in “TESTNEW” mode, which selects the best-fit evolution mode and conducts the phylogenetic analysis separately for each data set (Trifinopoulos et al. 2016). The CAT-GTR site mixture model implemented in PhyloBayes-MPI 1.7a (Lartillot et al. 2007) allows for site-specific rates of mutation, which is considered to be a more realistic model of amino acid evolution, especially for large multigene alignments (Maddock et al. 2016). PhyloBayes was run on the beta version of the Cipres server (https://cushion3.sdsc.edu/portal2/tools.action) (Miller et al. 2010), with default parameters (burnin = 500, invariable sites automatically removed from the alignment, two MCMC chains), and the analysis was stopped when the conditions considered to indicate a good run were reached (PhyloBayes manual: maxdiff <0.1 and minimum effective size >300). GHOST model is an edge-unlinked mixture model consisting of several site classes with separate sets of model parameters and edge lengths on the same tree topology, thus naturally accounting for heterotachous evolution (Crotty et al. 2017); analyses were run with 50,000 ultrafast bootstrap replicates (Minh et al. 2013). MLGO (Hu et al. 2014) was used to reconstruct the phylogeny using the gene-order (GO) data, with 1,000 bootstrap replicates, and an input file generated by PhyloSuite. Phylograms and gene orders were visualized in iTOL (Letunic and Bork 2007), and annotated using files generated by PhyloSuite. Skews were calculated and plotted using PhyloSuite and GraphDNA (Thomas et al. 2007).

Results

Models, Partitioning, and Compositional Heterogeneity Tests

Best-fit model for the nonpartitioned NUC data set was GTR+I+G4, mtZOA+F + I+G4 for the AAs, and SYM+R4 for the 18S. Best inferred partitioning strategy differed between NUC (12 partitions) and AAs (6 partitions) data sets (see supplementary file S1, Supplementary Material online, for details). The χ2 compositional homogeneity test denotes a sequence as “failed” if its nucleotide composition significantly deviates from the average composition of the alignment; in the nonpartitioned NUC data set only Penaeus vannamei (outgroup) passed the test (30 sequences failed); in the AAs data set 11 species failed; and in the 18S data set only Gammarus troglophilus (outgroup) failed.

NUC Data Set

GHOST, BI, and ML analyses (both partitioned and nonpartitioned) of the NUC data set produced highly congruent topologies (referred to as the “NUC-consensus” topology henceforth). Statistical support values were very high in BI (fig. 1; all inferred topologies available in the supplementary file S3, Supplementary Material online), relatively high in GHOST, and intermediate in ML analyses (low to high). Decapoda were rendered paraphyletic by the Euphausiacea nested within the clade (not in the ML-partitioned tree), and Mysida+Amphipoda were resolved as the sister-clade to Isopoda (see table 2 for the overview of some key features). In the isopod clade, Cymothoidae+Corallanidae families (Cymothoida) formed the basal clade, followed by Asellota and Phreatoicidea branches. The remaining isopods were divided into two sister-clades: Oniscidea and a “catch-all” clade comprising Limnoriidea, Valvifera, Sphaeromatidea, L. oceanica (a “rogue” Oniscidea species), and the remaining three Cymothoida species (Gyge ovalis, Bathynomus sp., and E.pulchra), which did not cluster together (fig. 1). Parsimony analysis produced a slightly rearranged topology, but crucial features were identical (table 2). Notable differences were: Asellota forming a sister-clade with G. ovalis + Lim. quadripunctata, and Phreatoicidea at the base of the catch-all clade. The PB analysis produced a notably different overall topology, with all nonisopod lineages forming a sister-clade to the isopods (monophyletic Decapoda), but the topology of the isopod clade was relatively similar to the NUC-consensus, apart from the rogue G.ovalis (fig. 2).

. 1.

—Mitochondrial phylogenomics of Isopoda (suborder information shown) reconstructed using partitioned nucleotide sequences of PCGs and rRNAs (NUC data set) and BI algorithm. A set of nine nonisopod Malacostraca species and Limulus polyphemus were used as outgroups (order information shown). The scale bar corresponds to the estimated number of substitutions per site. Bayesian posterior support values are shown next to corresponding nodes. Star sign indicates a putative origin of replication inversion scenario implied by the topology (see Discussion).

Table 2

Summary of Some Key Features of Phylogenetic Analyses

NUC	BI_P	BI_NP	ML_P	ML_NP	GHOST	Parsimony	PB_CAT
Basal	C+C	C+C	C+C	C+C	C+C	C+C	C+C
Cymothoida	Par	Par	Par	Par	Par	Par	Par
Sister-gr.	M+A	M+A	M+A	M+A	M+A	M+A	Malacostraca
Cym+Asel	+	+	+	+	+	+	+
Limnoriidea	G.o.	G.o.	G.o.	G.o.	G.o.	G.o.	*
AAs	BI_P	BI_NP	ML_P	ML_NP	GHOST	Parsimony	PB_CAT
Basal	Phreat.	C+C	Phreat.	Phreat.	Phreat.	C+C	Phreat.
Cymothoida	Par	Par	Par	Par	Par	Par	Par
Sister-gr.	M+A	M+A	Mysida	Mysida	Mysida	M+A	Mysida
Cym+Asel	+(s)	+(s)	+(s)	+(s)	+(s)	+	−
Limnoriidea	G.o.	G.o.	G.o.	G.o.	G.o.	G.o.	G.o.*
18S	BI		ML		GHOST	Parsimony	PB_CAT
Basal	Asellota		Asellota		Asellota	Limnoriidea	Asellota
Cymothoida	Par		Par		Par	Mon	Par
Sister-gr.	Decapodaa		Decapodaa		Decapodaa	Amphipoda	Amphipoda
Cym+Asel	−		−		−	−	−
Limnoriidea	Amphipodaa		Amphipodaa		Corallanidae	Basal	Cymothoida

Note.—Basal, basal isopod clade; Cymothoida, monophyletic or paraphyletic; Sister-gr., sister-group to the Isopoda; Cym+Asel, does the topology exhibit the artifact of a close relationship of Cymothoidae+Corallanidae and Asellota (isopods with congruent skews), where +(s) indicates a sister-group relationship; Limnoriidea, sister-group to this taxon, where G.o. is Gyge ovalis and * indicates polyphyly; P, partitioned; NP, nonpartitioned data set; C + C, Cymothoidae+Corallanidae; M + A, sister-group Mysida+Amphipoda; Phreat., Phreatoicidea.

Amphipoda clustered within the Isopoda clade (paraphyletic isopods).

. 2.

—A phylogram reconstructed using nonpartitioned NUC data set and an algorithm designed to address compositional heterogeneity: CAT-GTR (PB). Posterior Bayesian support values are shown. See figure 1 for other details.

Single-Gene Mitochondrial Data Sets

As we hypothesized that the topological instability may be driven by conflicting signals produced by different genes (Meiklejohn et al. 2014), we conducted ML phylogenetic analyses on 15 single-gene data sets (13 PCGs + 2 rRNAs). All these produced unique topologies (supplementary file S3, Supplementary Material online). Atp6 resolved Asellota+Cymothoidae+Corallanidae as the basal isopod clade; atp8 (a very small gene) produced an almost nonsensical topology (defined as: in stark disagreement with any reasonable phylogenetic hypothesis); nad1 resolved Asellota at the base; nad2 produced a slightly rearranged NUC-consensus topology; nad3 (small gene) a nonsensical topology; nad4 produced a rearranged topology (compared with NUC-consensus); nad4L (small gene) produced a highly rearranged isopod clade; despite its large size (≈1,700 bp), nad5 produced a noncanonical isopod topology; nad6 produced a nonsensical topology (paraphyletic Isopoda, rogue G. ovalis); cox1 produced a highly rearranged topology, including both the nonisopod malacostraca and isopod clades, with Phreatoicidea+Cymothoidae+Corralanidae as the basal isopod clade, and noncanonical Valvifera position; cox2 produced a slightly rearranged NUC-consensus topology; cox3 a highly rearranged topology, with a unique Asellota+Corallanidae clade at the isopod base; cytb produced basal Phreatoicidea, and the remaining taxa divided into Oniscidea and “catch all” sister-clades, where the rogue Cymothoidae+Corallanidae (along with Asellota) were on a long branch in the derived part of the clade; rrnS (or 12S) topology was in some aspects similar to the cytb, but with Phreatoicidea+Limnoriidea as the basal isopod clade; rrnL (or 16S) produced an almost nonsensical topology, with paraphyletic Isopoda and rogue G. ovalis.

Amino Acids Data Set

Amino acids (concatenated 13 PCGs) produced topologies that differed from the NUC-consensus one, with instable topology of the nonisopod Malacostraca, including paraphyletic Eucarida (Decapoda+Euphausiacea) in five (BI/ML, nonpartitioned/partitioned, and GHOST; fig. 3) out of six (PB; fig. 4) analyses. Sister-clade to Isopoda also varied: Mysida+Amphipoda or Mysida (table 2). In the isopod clade, partitioning had a major effect on the BI analysis: nonpartitioned data set resolved Cymothoidae+Corallanidae at the base, followed by Phreatoicidea+Asellota. Parsimony analysis produced a similar topology, but Phreatoicidea and G. ovalis + Lim. quadripunctata clades switched places. The remaining five analyses (including the partitioned BI data set) resolved Phreatoicidea as the basal clade (table 2), and the sister-group to the remaining isopods (minus Phreatoicidea) was Asellota+Cymothoidae+Corallanidae in four analyses, and Asellota in PB. Five analyses (minus PB) produced a topology of the remainder of the isopod clade that was partially congruent with the NUC-consensus topology, but instead of it being divided into Oniscidea + all other taxa, here G. ovalis+Lim. quadripunctata were at the base (except in Parsimony: Phreatoicidea). Oniscidea (rendered paraphyletic by L. oceanica) clade topology was stable in all six, but the topology of the catch-all clade was different from the NUC-consensus topology.

. 3.

—A phylogram reconstructed using amino acid data set (AAs; 13 PCGs) in combination with data partitioning strategy and BI algorithm. See figure 1 for other details.

. 4.

—A phylogram reconstructed using AAs data set and CAT-GTR algorithm designed for heterogeneous data sets (PB). See figure 1 for other details.

PB produced the only topology (fig. 4) with monophyletic Oniscidea, with the rogue L. oceanica at the base of the clade. Also importantly, Cymothoidae+Corallanidae clade was placed on a long branch within the catch-all clade, together with E.pulchra (all Cymothoida). The suborder was still rendered paraphyletic by the positions of Bathynomus sp. and rogue (polytomy) G. ovalis.

Gene Order-Based Phylogeny

As gene orders are unlikely to exhibit homoplasy, it has been hypothesized that they may be able to resolve difficult (deep) phylogenies in some cases (Boore 2006), so we tested this approach. As many tRNA genes were missing or we suspected that they may be misannotated, to test for the presence of false signals, we used two data sets: one excluding all tRNAs (PCGs+rRNAs) and one (PCGs+rRNAs+tRNAs) containing all (identified) tRNAs apart from the ambiguously annotated trnL genes. Both data sets produced a number of paraphyletic major clades and largely nonsensical topologies with very low bootstrap support values (supplementary file S3, Supplementary Material online: GO trees). Furthermore, the two data sets produced incongruent topologies, and different runs of the PCGs+rRNAs data set produced different topologies (supplementary file S3, Supplementary Material online).

Nuclear Marker-Based Phylogeny (18S)

To view the phylogeny of Isopoda from a nonmitochondrial perspective, we used the nuclear 18S gene. This approach can also help us test whether the underlying reason for the conflicting phylogenetic signals between different studies might be mitochondrial introgression, which, recent evidence shows, is more widespread than previously thought (Jakovlić et al. 2013; Edwards et al. 2016; Mallet et al. 2016). Results were also unstable, that is, different runs and data sets (addition and removal of taxa) would often produce different topologies, but most of them produced Asellota at the base (Limnoriidea in Parsimony, but with very low support); highly derived Cymothoida, rendered paraphyletic by the nested Limnoriidea (monophyletic in Parsimony); and Oniscidea rendered paraphyletic by the rogue Ligiidae clade (Ligia and Ligidium genera). BI, ML, and GHOST analyses produced relatively congruent topologies, with Isopoda rendered paraphyletic by the Amphipoda nested within the Cymothoida (supplementary file S3, Supplementary Material online). As PB and Parsimony analyses produced monophyletic Isopoda, with Amphipoda as the sister-clade (fig. 5), we can confidently reject this as a compositional heterogeneity artifact. Cymothoida were monophyletic only in the Parsimony analysis, and divided into two clades in all topologies: a stable monophyletic clade comprising Dajidae and Bopyridae; and a large instable (exhibiting pervasive paraphyly) clade comprising Cirolanidae, Corallanidae, and Cymothoidae families, and aforementioned “intruders” (Limnoriidea and Amphipoda). Corallanidae and Cirolanidae were mostly paraphyletic (somewhat erratic behavior of E.pulchra), whereas Cymothoidae (monophyletic) were highly derived and exhibited disproportionately long branches.

. 5.

—A phylogram inferred using the nuclear 18S gene and CAT-GTR algorithm (PB). See figure 1 for other details.

Discussion

Apart from the heterotachous GHOST model, which tended to produce results identical to the common ML analysis, all other variables produced notable impacts on the topology (table 2), with unique topologies by far outnumbering identical topologies. As regards the major unresolved issue of the isopod phylogeny, the sister-clade to all other isopods (basal clade) and the monophyly of Cymothoida, mitochondrial nucleotides (NUC) consistently produced Cymothoidae and Corallanidae as the basal clade, first followed by Asellota, then by Phreatoicidea, whereas other Cymothoida (paraphyletic) tended to be scattered throughout the central “catch-all” clade. AAs, however, largely resolved Phreatoicidea as the basal isopod clade, with all remaining isopods split into two sister-clades: 1) Asellota + Cymothoidae/Corallanidae and 2) all remaining taxa. The Parsimony method produced Cymothoidae+Corallanidae at the base, followed by Asellota + (G. ovalis + Lim. quadripunctata) using both data sets. PB analysis of AAs data set produced a remarkably different topology, with Phreatoicidea at the base, followed by Asellota, but Cymothoidae+Corallanidae were relatively derived, and Cymothoida not so deeply paraphyletic. Nuclear (18S gene) topology was also strongly affected by the methodology, but mostly resolved Asellota as the basal isopod clade (Limnoriidea in Parsimony), and Cymothoida as paraphyletic, with nested Limnoriidea (not in Parsimony), but highly derived. Although it appears that we did not manage to reach a conclusion, we did manage to identify a feature of isopod mitogenomes that may explain this instability.

Base Composition Skews

Organellar genomes often exhibit a phenomenon known as strand asymmetry, or strand compositional bias, where positive AT skew values indicate more A than T on the strand, positive GC skews indicate more G than C, and vice versa (Reyes et al. 1998; Wei et al. 2010). This is believed to be caused by hydrolytic deamination of bases on the leading strand when it is single stranded, that is, during replication and/or transcription (Reyes et al. 1998; Bernt et al. 2013; Fonseca et al. 2014). Whereas other crustacean taxa usually exhibit positive overall AT skews for genes located on the plus (majority) strand and negative GC skews for genes on the minus (minority) strand (Hassanin 2006; Wei et al. 2010), isopod mitogenomes usually exhibit an inverted skew pattern (Kilpert and Podsiadlowski 2006; Kilpert et al. 2012; Yu et al. 2018). This is believed to be a consequence of an inversion of the replication origin (RO), where the changed replication order of two mitochondrial DNA strands consequently resulted in an inversed strand asymmetry (Hassanin et al. 2005; Kilpert and Podsiadlowski 2006; Wei et al. 2010; Kilpert et al. 2012; Bernt et al. 2013).

It is known that Asellus aquaticus (Asellota) possesses an inversed skew in comparison to other isopod taxa (Kilpert and Podsiadlowski 2006), but here we found that the three available Cymothoidae and Corallanidae species (C. indica, T. chinensis, and I. japonensis) also exhibit inversed skew patterns (fig. 6 and table 1). As regards nonisopod mitogenomes in the data set, they exhibit GC skews comparable to the isopod outliers, from −0.014 in Metacrangonyx repens to −0.332 in Typhlatya miravetensis, and mixed AT skews (table 1). The inversed skew in A.aquaticus led Kilpert and Podsiadlowski (2006) to speculate that this taxon branched off first in the isopod phylogeny, but a few years later Kilpert et al. (2012) noticed that the basal position of Phreatoicidea causes a conflict in explaining the inversed isopod GC skew (present in Eophreatoicus sp.; fig. 4—red stars). Our NUC data set would mostly support a modified version of the first scenario, with the RO inversion occurring in isopods after Cymothoidae+Corallanidae and Asellota branched off (figs. 1 and 2—red stars), whereas AAs would support an RO inversion in the common ancestor of Asellota+(Cymothoidae+Corallanidae) clade (fig. 3). Parsimony analyses further complicate this by Asellota forming a sister-clade with G. ovalis and Lim. quadripunctata. However, there is no support for either of these scenarios from nuclear (18S) or morphological (Wilson 2009) data, which relatively consistently indicate that Asellota (18S) or Asellota+Phreatoicidea (morphology) are the basal clade.

. 6.

—Cumulative GC skews of the majority strands of a selected subset of mitogenomes used for phylogenetic analyses.

Variations in base composition can bias phylogenetic analyses (Romiguier and Roux 2017), and skew-driven LBA phylogenetic artifacts have been reported in arthropods (Hassanin et al. 2005; Hassanin 2006) and other metazoans (Sun et al. 2018). Most mitochondrial analyses produced topologies where taxa with congruent skews clustered together: outgroups+(Cymothoidae+Corallanidae)+Asellota (negative GC skews taxa), followed by the remaining Isopod taxa (positive GC skews). Therefore, we can conclude with high confidence that mitochondrial data produce artifactual clustering of branches exhibiting similar skews: the inversed skew of highly derived Cymothoidae+Corallanidae results in their clustering at the base of the isopod clade, phylogenetically close to other taxa with similar (homoplastic) skew patterns: Asellota and nonisopod Malacostraca. This explains the topological instability, incongruent phylogenetic hypotheses, and it rejects the above RO inversion scenarios. We can conclude that inversed skews of Asellota and Cymothoidae+Corallanidae are nonsynapomorphic, and that highly derived Cymothoidae+Corallanidae underwent an additional RO inversion, which eventually resulted in a double-inverted skew (homoplastic with Asellota and other Malacostraca).

In support of this, the PB (heterogeneous) analysis of the AAs data set produced an mtDNA topology that exhibited notable similarity to the 18S and morphology-based topologies, where Cymothoidae+Corallanidae clustered with E. pulchra (Cirolanidae) in the relatively derived part of the isopod clade. This indicates that a combination of the AAs data set, which is expected to be less affected by skews than nucleotides, and a heterogeneous CAT-GTR model, was the most successful in attenuating the phylogenetic artifacts caused by compositional biases. It also confutes the hypothesis that conforming skews between Asellota and Cymothoidae+Corallanidae may be a result of a mitochondrial introgression event.

Collation of Evolutionary Hypotheses

Having established the fact that taxa with conforming homoplastic skews cluster together in both AAs and NUC mitochondrial data sets (table 2), we can infer the most parsimonious hypothesis for the course of events in the evolutionary history of Isopoda. First, we can reject with confidence the basal position of Cymothoidae+Corallanidae as an artifact. This indicates that the basal isopod taxon is either Asellota (Wilson 2009), Phreatoicidea (Brusca and Wilson 1991), or Asellota+Phreatoicidea sister-clade (Wilson 1999; Dreyer and Wägele 2001; Kilpert et al. 2012). The latter two scenarios are less parsimonious, as they would require at least three independent RO inversions in the evolutionary history of Isopoda (in the ancestral isopod, in Asellota, and in Cymothoidae+Corallanidae), whereas the first scenario is more parsimonious, as it requires only two (in the ancestral isopod after the split of Asellota and in Cymothoidae+Corallanidae; fig. 5—red stars). Additionally, the 18S data set relatively consistently resolved Asellota as the basal branch (disregarding paraphyletic Isopoda and Parsimony analysis; table 2). Therefore, we can tentatively conclude that multiple evidence supports the original hypothesis of Kilpert and Podsiadlowski (2006): Asellota is the oldest isopod branch and RO inversion in isopods occurred after the Asellota branched off.

Although this resolves the issue of the deep paraphyly of Cymothoida, that is, places the rogue Cymothoidae+Corallanidae clade back within the remaining Cymothoida, 18S data still resolve the Cymothoida paraphyletic, with Limnoriidea nested between two large clades (Dajidae+Bopyridae and Cirolanidae+Corallanidae+ Cymothoidae). As Limnoriidea was resolved as the basal isopod taxon in the 18S Parsimony analysis, and Corallanidae+Cymothoidae exhibit a disproportionately long branch, we suspect that this is an LBA between two taxa exhibiting elevated evolutionary rates. The monophyly of Corallanidae is unsupported by our 18S analyses, so it will be necessary to sequence further mitogenomic (to identify skews) and nuclear data for this group of cymothoid families. This combination of skews and nuclear data would enable us to identify the exact point in the evolutionary history where the RO inversion occurred in these taxa, and infer the most parsimonious topology and/or taxonomy, that is, the one that supports a single RO inversion (or introgression event), as opposed to those that would require multiple events.

Our analyses further corroborated the existence of several rogue taxa that exhibit somewhat erratic topological behavior. The position of L.oceanica (nominally Oniscidea: Ligiidae), a recognized rogue taxon (Wilson 2009; Lins et al. 2017; Shen et al. 2017; Yu et al. 2018), was mostly resolved at the base of the “catch-all” clade in mtDNA analyses. However, in the (putatively) most reliable mitochondrial topology AAs+PB, it was resolved as the basal Oniscidea species. Although this is in perfect agreement with morphological data, which resolve Ligiidae as the most primitive Oniscidea clade (Schmidt 2008), we cannot claim that this issue is fully resolved, since the entire Ligiidae family exhibited rogue behavior in the 18S data set as well. Three Cymothoida taxa in the mtDNA data set that exhibit standard isopod skews, Bathynomus sp., E. pulchra (both Cirolanidae), and G. ovalis (Bopyridae), also exhibited rather instable topological behavior, and we did not find support for their monophyly using the mtDNA data. Although the Cirolanidae are believed to be ancient (Wetzer 2002) and highly plesiomorphic within this suborder (Brandt and Poore 2003), this is not supported by the 18S data set. As all three species exhibit highly rearranged gene orders (supplementary file S3, Supplementary Material online: Gene arrangements), and as there is evidence of a close positive correlation between the mitogenomic architectural instability and the mutation rate (Shao et al. 2003; Hassanin 2006; Xu et al. 2006), we hypothesize that frequent genome rearrangements may have resulted in an accelerated mutational rate in these species. In agreement with this hypothesis, another rogue species, Lim. quadripunctata (Limnoriidea) (Wilson 2009; Lloyd et al. 2015; Lins et al. 2017; Zou et al. 2018), also exhibits a highly rearranged gene order (supplementary file S3, Supplementary Material online: Gene arrangements), and E. pulchra and Lim. quadripunctata exhibit some of the highest evolutionary rates among the isopod (and decapod) mitogenomes (Shen et al. 2017). There are other indications that Limnoriidea is evolving under unique evolutionary pressures: Lim. quadripunctata clustered with G. ovalis in most of our mitochondrial analyses (table 2) and it was resolved as the basal isopod clade in a previous mitogenomic study (Cymothoidae+Corallanidae unavailable at the time) (Lins et al. 2017). Limnoriidea was the basal isopod clade in our 18S Parsimony analysis, and it nested within the Cymothoida in the remaining 18S analyses (sometimes together with Decapoda) (table 2). Furthermore, inverted mitogenomic skews also do not explain the disproportionately long branch of Cymothoidae in the 18S data set. This is an indication that there are other factors (aside from the RO inversions) that cause accelerated substitution rates in Cymothoidae, Corallanidae, and Limnoriidae, so we urge further studies of these taxa, as they appear to be exceptionally interesting from the molecular evolution perspective.

The topology of nonisopod Malacostraca was rather instable as well, with PB producing notably different topologies from other analyses. As nonisopod Malacostraca also exhibit notable variability in skews (table 1), we conclude that compositional heterogeneity also interfered with phylogenetic reconstruction. The relevant question for this study is that of the sister-group to Isopoda: Mysida+Amphipoda in most of our mtNUC analyses (except PB), BI+AAs, and all Parsimony topologies; Mysida in most AAs analyses (ML+GHOST+PB); and Amphipoda in the PB 18S analysis (Amphipoda were nested within the Isopoda in ML and BI 18S topologies) (table 2). As Amphipoda are considered to be the most prominent contender for this position (Wilson 2009), we can conclude that our analyses support this close relationship.

Methodological Implications

Our findings show that mitochondrial sequence data are producing artifactual LBA relationships in isopods, and thus are a poor tool for the reconstruction of their phylogeny. Intriguingly, these did not affect only the nucleotide data set, but also the amino acid data set, which should be less affected by nonadaptive compositional biases, as nonsynonymous mutations are likely to be affected by the purifying selection (as opposed to synonymous mutations). This corroborates that mitochondrial strand asymmetry (skews) can have very pronounced effects on the composition of encoded proteins (Min and Hickey 2007; Botero-Castro et al. 2018). It should be noted that a more thorough sampling of mitogenomes may have resulted in somewhat different artifacts; for example, had we included only one outgroup and one Asellota mitogenome, and a large number of Cymothoidae+Corallanidae, it is likely that we would have obtained a topology with the derived C + C, and Asellota (and possibly even the outgroup as well), clustering within the Cymothoida (inversed direction of the skew “gravity”). Regardless of the data set, it remains highly likely that attraction between the homoplastic skews of these two (distant) clades would confound the phylogenetic analysis.

All individual genes produced unique topologies, which has important implications for the interpretation of previous results inferred using single-gene data sets. This phenomenon has been observed in isopods before on a much smaller scale (Wetzer 2002), and it is in agreement with the proposed mosaic nature of (mitochondrial) genomes (Pollard et al. 2006; Degnan and Rosenberg 2009; Romiguier et al. 2013; Romiguier and Roux 2017). Intriguingly, although 12S gene produced a topology that placed Cymothoidae+Corallanidae in the derived part of the clade, Asellota clustered within the Cymothoida clade (supplementary file S3, Supplementary Material online), which supports the observation that taxa with conforming skews tend to cluster together regardless of the direction of artifactual “gravity-pull.” As gene orders produced topological instability, very low support, and almost nonsensical topologies, with highly rearranged orders at the base of the isopod clade, we hypothesize that discontinuous evolution of mitogenomic architecture evolution (Zou et al. 2017) produces phylogenetic artifacts such as LBA, and renders them useless for the task.

All of the tested standard models (BI, ML, and Parsimony) were very sensitive to compositional biases, and produced strong artifacts. This included the new heterotachous model, GHOST, which mostly produced results identical to the ML algorithm. Importantly, as regards the aforementioned feud about the most suitable methodological approach to account for compositional heterogeneity (Feuda et al. 2017; Whelan and Halanych 2017), our results indicate that (in isopods) the PB CAT-GTR model by far outperforms the partitioning (assigning different evolutionary models to different partitions). Although we discourage the use of mitochondrial data as a phylogenetic tool for the reconstruction in isopods (even in combination with nuclear and morphological data), they may still be of some use in studies focused only on the isopod suborders that exhibit congruent skews. It should be mentioned that there are other available methodological approaches designed to account for compositional heterogeneity (p4, PHASE, nhPhyML and so on) (Foster 2004; Boussau and Gouy 2006; Hassanin 2006; Gowri-Shankar and Rattray 2007; Sheffield et al. 2009; Richards et al. 2018; Yang et al. 2018), so future studies may attempt to test their performance as well.

Conclusions

With respect to our working hypothesis, we can accept the first part of it: asymmetrical mutational pressures do generate compositional heterogeneity in isopod mitogenomes and interfere with phylogenetic reconstruction. However, we were mistaken in assuming that these pressures are primarily adaptive, that is, caused by their radically diverse life histories. Instead, they appear to be driven largely by the origin of replication inversions. These appear to generate very strong compositional biases, which render the mitogenomic sequence data a very poor tool for reconstructing the evolutionary history of Isopoda. None of the tools used here managed to fully revolve these compositional biases, but PB CAT-GTR algorithm outperformed partitioning, and best results were achieved by combining it with the amino acids data set. As mtDNA data have played a major role in our understanding of the evolutionary history of life on Earth (Rubinoff et al. 2005), implications of this study are much broader than its original scope. As our findings indicate that architectural rearrangements might produce major compositional biases even on comparatively short evolutionary timescales, proving the suitability of data via GC and AT skew analyses should be a prerequisite for any study that aims to use mitochondrial data for phylogenetic reconstruction, even among the closely related taxa. These findings should not discourage scientists from sequencing further isopod mitogenomes (especially interesting are nonrepresented Asellota and Cymothoida families), as their architectural hypervariability still makes them a useful tool for unraveling the conundrums of evolution of mitochondrial architecture, and as mitochondrial skews can be used as an additional phylogenetic tool to infer the most parsimonious phylogenetic hypotheses.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.