New insights into Oniscidea (Crustacea: Isopoda) mitogenome structural features and phylogenetic placement of targeted taxa using mitogenomic and nuclear data

Among Isopoda suborders the Oniscidea has the highest species richness, and is also the largest terrestrial group in the Crustacea. Terrestrial isopods are an excellent case to study adaptations related to sea-land transition. However, the monophyly of Oniscidea and the relationships of the main lineages has been debated over the last three decades. Aiming to further explore structural features of mitochondrial genome and investigate the phylogenetic relationships within Oniscidea, the mitogenomes and a series of nuclear markers of the oniscids Ligia exotica and Mongoloniscus sinensis were sequenced. The nuclear genome was represented by four nuclear genes analyzed in a separate dataset. The mitogenomes of L. exotica and M. sinensis were 16,018 and 14,978 bp in length, respectively. Both included 13 protein-coding genes, 2 rRNA genes, and 21 and 19 tRNA genes respectively, missing one and three tRNA genes respectively compared to the isopod ground pattern. The M. sinensis mitogenome had higher average A+T content (~75.3%) than any other isopod studied to date. Phylogenetic analyses conrmed the assignment of M. sinesis to Agnaridae, as well as the sister-group relationship of the family with Porcellionidae, one of the more derived Crinochaeta clades. On the other hand, the basal position of Ligia within Oniscidea and the close evolutionary relationship with the aquatic groups Valvifera, Shaperomatida and some Cymothoidea that were included in our analysis, is indicated. Aiming to further explore oniscid mitogenome characters and investigate the phylogenetic relationships within the group and between terrestrial isopods and aquatic relatives, the complete mitogenomes of L. exotica and M. sinensis were sequenced. Furthermore, we also sequenced 18s, 28s ribosomal RNA genes, protein-coding Sodium-Potassium Pump (NAK) and Phosphoenolpyruvate Carboxykinase (PEPCK) genes of M. sinensis and Ligia cinerascens to supplement the assignment based on mitogenome dataset. Mitogenome and nuclear data were kept separately taking into account the limited overlap in the taxa and the fact that they originate from organelles under different evolutionary pressure.


Introduction
The Isopoda, with more than 10,300 described species, is the largest order within the Peracarida, and includes an amazing diversity in morphology and ecology [1,2]. Isopods are diverse in the ocean, land, and freshwater, include free-living to parasitic taxa, and range from the deep oceans to at least 4,700m above sea level. Currently, the Isopoda is divided into 11 suborders, with Oniscidea being the most speciose suborder consisted of 39 families, 523 genera and more than 3700 species [3][4][5]. From an ecological point of view, although Oniscidea can be found across most terrestrial biomes, from nearshore marine habitats to moist forests and arid deserts, they exhibit a strong preference for humid micro-habitats [5]. These phytosaprophagous animals are among the most important decomposers [6]. The suborder is currently separated into ve lineages: Diplocheta, Tylida, Microcheta, Synocheta and Crinocheta [7][8][9]. The most basal members of the group, in the families Ligiidae, Ligidiidae, and Tylidae are amphibious and restricted to littoral marine habitats [10]. Ligia Fabricius 1798 is of particular interest as it is considered to be an intermediate form with morphological, physiological and behavioral characteristics that exhibit similarities with marine ancestors and fully terrestrial Oniscidea [8,11]. Schmidt (2008) in a detailed morpho-phylogenetic study of oniscids argued that the Ligiidae was the most basal clade in the Oniscidea, suggesting a marine origin for the suborder. Hence the phylogenetic relationship of Ligiidae with the rest of Oniscidea and aquatic relatives is essential in our effort to understand the evolutionary steps towards Isopods terrestrialization.
Extremely high genetic distances in both nuclear and mitochondrial genome reaching up to 50.3% between confamiliar genera and 20.3% between individuals of the same species were identi ed [24,25]. These results highlight the vivid divergence between even closely related taxa indicating the need of phylogenetic assessment using genetic data where possible.

Nuclear DNA primer design, LA-PCR ampli cation and sequencing
Primers targeting the 18s rRNA, 28s rRNA, Potassium Pump (NAK) and Phosphoenolpyruvate Carboxykinase (PEPCK) genes were designed according available sequences of closely related species (S1 Table). PCR conditions were as follows: initial denaturation 94°C for 2 min, then 35 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 1 min/kb, followed by the nal extension at 72°C for 10 min. The total volume for PCR and LA-PCR was 50 μl, of which Takara LATaq (5 U/μl) was 0.5 μl, 10×LATaq Buffer II (Mg2+) was 5 μl, dNTP mixture (2.5 mM) was 8 μl, template was 60 ng, and the total volume was reached with distilled water. The nal concentration of the forward and reverse primers was 0.2~1.0 μM, and that of MgCl2 was 2.0 mM. The PCR products were sequenced directly, or if needed rst cloned into a pMD18-T vector (Takara, JAP) and then sequenced, by the dideoxynucleotide procedure, using an ABI 3730 automatic sequencer (Sanger sequencing) using the same set of primers.

Phylogenetic analysis
Phylogenetic analyses were conducted using concatenated amino acid sequences of all 13 protein-coding genes (PCGs) from 2 new sequenced and 28 currently available isopod mitogenomes from GenBank. Amino acid sequences were aligned using MAFFT in batch mode under L-INS-i strategy [44,45]. Aligned loci were concatenated with PhyloSuite. Bayesian analysis using the probabilistic CAT-GTR model was implemented in PhyloBayes-MPI 1.7a through CIPRES server [46]. The analysis run using default parameters until maxdiff < 0.1 and minimum effective size > 300 were reached. Multiple sequence alignments were performed using MAFFT wih Q-INS-i algorithm for not-coding 18s rRNA and 28s rRNA and '--auto' option for NAK and PEPCK. Ambiguously or poorly aligned regions were removed using Gblocks allowing smaller nal blocks, less strict anking positions and presence of gaps at the nal alignment [47,48]. The best-t partitioning scheme and evolutionary model for each partition were selected with PartitionFinder2 under BIC criterion [49]. The nal concatenated dataset set was partitioned by gene into four distinct data blocks. Bayesian inference analysis (BI) was implemented using MrBayes v.3.2.6 [50], with the selected model of nucleotide evolution for each gene, allowing rate heterogeneity between partitions.

Complete mitochondrial genomes
The complete mitochondrial genome of M. sinensis and L. exotica are circular, double-stranded DNA molecules, 16,018 and 14,978 bp in length, respectively. The M. sinensis mitogenome is comprised of 13 protein-coding genes, 2 rRNA genes and 19 tRNAs, with trnA, trnE and trnL1 are missing. In addition, 18 overlapping regions and 13 intergenic regions were identi ed ( Table 2). The average A+T content of the M. sinensis mt genome is approximately 75.3% (S2 Table), higher than any other known isopod mitogenomes (range: 54.4 -71.2 ). The mitogenome of L. exotica is composed of 13 PCGs, 2 rRNA genes and 21 tRNA genes, only lacking the trnG gene. In total 17 overlapping regions and 14 intergenic regions were detected ( Table 3). The average A+T content is approximately 59.1% (S2 Table). The PCG size and start/stop codons of both species is typical of other isopod mitogenomes (Tables 2 and 3). Consistent with the inversed strand asymmetry shown in most isopods, GC skews of the mitochondrial genomes of 12 examined Oniscidea was positive (S2 Table).

Transfer RNA genes
Nineteen tRNAs were found in M. sinensis, ranging from 51 bp (trnS2) to 73 bp in size (trnS1), totaling 1158 bp. In L. exotica, 21 tRNAs were identi ed, ranging from 53 bp (trnC) to 64 bp in size (trnI, trnM, trnW), summing up to 1278 bp. Transfer RNA genes are distributed throughout the mitogenome and found on both strands. The proposed secondary structures of all identi ed tRNAs are shown in Figs 1 and 2. The majority of tRNAs have a common t-shaped or cloverleaf secondary structure. In M. sinensis exceptions include the absence of DHU-arm in trnC, and TΨC-arm in trnG, trnK, trnP, trnS2, trnW and trnV. In L. exotica exceptions are the absence of DHU-arm in trnC and trnV, and TΨC-arm in trnF, trnP, trnM, trnL2, trnK, trnD, trnR.

Gene translocations of mitogenomes
The comparison of mitogenomes of M. sinensis and L. exotica with the known ground pattern of isopods revealed four and three gene rearrangements respectively (Fig 3). In M. sinensis i) trnR is located between cox3 and nad3 instead of trnA and nad1 where it is usually found; ii) trnV is between nad3 and nad1, instead of between 16s rRNA and trnQ, iii) trnG is between trnW and trnT, instead of between cox3 and nad3 and iv) trnT is located between trnG and trnS1, instead of between cob and nad5. Similarly, in L. exotica, trnR was translocated between cox3 and nad3, while trnW and trnE are interchanged. Compared with the known found pattern of isopods, both M. sinensis and L. exotica have two and one gene missing respectively. In M. sinensis i) trnA is missing instead of between nad3 and trnR; ii) trnL1 is missing instead of between nad1 and trnN. In L. exotica, trnG is missing instead of between cox3 and nad3.Gene translocations as well as loss of tRNA genes are well known in isopod mitogenomes. Some tRNA genes are missing in every oniscids studied to date. Some of the genetic rearrangements are attributed to the missing tRNA genes.

Amino Acids Dataset (13PCGs)
Phylogenetic reconstruction based on the concatenated amino acid sequences of all 13 PCGs recovered a statistically well supported monophyletic Oniscidea clade, with Ligiidae in basal position. M. sinensis was recovered as sister to a clade to Porcellionidae represented by Porcellionides and Porcellio (Fig 4).

Nuclear Dataset (18s-28s-NAK-PEPCK)
The aligned data for the combined nuclear genes was 1719 base pairs (bp). The best substitution models were SYM+I+G for 18s, GTR+G for 28s, SYM+I+G for NAK, and GTR+I+G for PEPCK. The Oniscidea was recovered in two separate statistically well supported clades, with Ligia species in one, and the rest of oniscids in another including Sphaeroma. Although the basal phylogenetic relationships of the group are not fully resolved, the close evolutionary relationship of Oniscidea excluding Ligia, with Sphaeroma is statistically supported. Typhloligidium Verhoeff, 1918, Ligidium Brandt, 1833 and Taurologidium Borutzky, 1950 all members of the newly established family Ligidiidae are grouped together [9], Agnaridae, Platyarthridae, Trichoniscidae, Trachelipodidae and Porcellionidae were all recovered as polyphyletic. Mongoloniscus sinensis was sister to Agnara madagascariensis and in turn with to three species of Hemilepistus, but in a separate branch than Protracheoniscus aff. fossuliger the only other agnarid included.

Characteristics of mitogenomes
The aim of this work was to determine the mitogenome structure of two Oniscidea species and exploit these data for the phylogenetic assignment of targeted species and possibly revisit phylogenetic relationships within Oniscidea and/or Isopoda. The mitogenomes of L. exotica and M. sinensis are circular double stranded, with incomplete or totally missing tRNA genes. A recently published work speculated that these tRNA gene fragments may actually be transcribed and edited into functional tRNAs, since their presence is common in isopods [23]. Structural rearrangements of these two species' mitogenomes may be explained by linearized and fragmented events [51].
Unfortunately, the pattern of speci c genes loss and rearrangement, and the exact mechanism behind it, had not been found in the 12 known species of the suborder Oniscidea, possibly due to insu cient data. It is hoped that the pattern, exact mechanism and the meaningful information behind the pattern can be found in future further studies.

Data Processing
Considering the confounding phylogenies resulting from different datasets based on nuclear or/and mitochondrial data [9,14,15,17,23], we compiled two different inclusive datasets in order to assess the patterns yielding from genetic loci with different evolutionary history. Herein we investigated the phylogenetic assignment of the targeted species within Oniscidea combining the results from separate analyses of genomic data from different organelles. Taking into account the long evolutionary history of the group, in order to avoid saturation events due to high mutation rates of mitogenomes, only PCGs were included in mitochondrial datasets analyses [52]. Furthermore, these loci were translated in AA sequences and run under a heterogeneous CAT-GTR model. In this way we avoided phylogenetic artifacts attributed to biases due to similar GC skews among distant taxa [22]. On the other hand, nuclear non-coding genes were treated with Gblocks in order to eliminate poorly aligned regions. With this treatment we excluded from our analyses DNA fragments accumulating mutations at higher rates, compared to their anking regions, which could result in misleading phylogenies.

Phylogenetic relationships of Oniscidea
Phylogenetic relationships within Oniscidea and other isopod groups, even at high taxonomic scales, based on morphological characters or molecular data were repeatedly debated [9,15,[53][54][55]. The monophyly of Oniscidea is well supported by morphological characters such as the complex water-conducting system and reduced rst antenna with only three articles [8,12,13,15,56,57]. Although it was alternatively suggested as a possible homoplasy, the rst antenna with rudimentary three articles has been regarded as a prominent synapomorphy for Oniscidea [8,12]. Furthermore, taking into account the existence of two types of water conducting systems (i.e "Ligia type" and "Porcellio type") [10], we could speculate that these two characters might not be attributed to the monophyletic origin of Oniscidea but to convergent evolution as a response to environment challenges [15].

Phylogeny based on mitogenome dataset
Monophyly of Oniscidea have been questioned by several phylogenetic analyses including molecular data [9,14,22,23,58]. Solely based on mitochondrial data recently published works failed to support the monophyly of Oniscidea [20,21,51]. The monophyly of Oniscidea mainly relies on the phylogenetic origin of the key genera Ligia, Ligidium, Tylos and Helleria whose close evolutionary relationship with marine ancestors was repeatedly highlighted in the past [9,15,22]. Therefore, two species of the genus Ligia were added in this paper to explore the monophyly of the suborder Oniscidea. Particular focus was given in Ligia, as that genus have been proposed to represent an intermediate form in the evolution of terrestriality in isopods [11]. However, the phylogenetic placement of Ligia within Oniscidea was debated in the past [8,9].
According to our mitogenome analysis results, Ligia appears to form a statistically well supported monophyletic group with Crinocheta, the more evolutionary recent Oniscidea lineage whose representatives exhibiting some of the most pioneer adaptations to terrestrial life. These results are also supported by a series of well described Oniscidea synapomorphies [7,8]. Based on the same mitochondrial loci, running analyses under different evolutionary substitution models recently published works came up with different results [15,23]. In contrast to Lins et al. (2017) our results are in agreement with Zou et al. (2020) ndings whose analysis was conducted under a CAT-GTR model (Fig 4).
Furthermore, our analysis revealed the close evolutionary relationship of Oniscidea with Sphaeromatidea, Valvifera and some Cymothoidea. Sphaeromatidae and Valvifera were also shown to be the more closely related marine taxa to terrestrial isopods by recently published works based on mitochondrial and nuclear data [9,23]. All Asellota members are grouped together with strong statistical support while Cymothoa although represented by a small number of species appear to be paraphyletic. Phreatoicidea should be considered as the more ancestral isopod lineage rooting the mitogenome tree.

Nuclear markers Phylogeny
Although the mitochondrial genome data were mainly used in this paper to elucidate the phylogenetic relationship of Oniscidea, nuclear gene data were also used for supplement, as detailed in the appendix materials. The main differences between the nuclear dataset presented by Dimitriou et al. (2019) and the one analyzed herein, mainly enriched with outgroups, are located at the basis of the constructed phylogenetic tree and added a key genus species, L. cinerascens.
The basal position of Phreatoicidea as well as Asellota at Isopods phylogeny is also supported by data of nuclear origin. On the other hand, we failed to save the monophyly of Oniscidea as Sphaeroma appear to be more closely related to terrestrial isopods than Ligia (S3 Fig.) (Hence the monophyly of Diplocheta could not be saved as well. Νumerus taxonomic revisions at lower taxonomic levels especially withing the more apical sister groups Crinocheta and Synochetawere proposed [5,24].
The familiar assignment of Mongoloniscus was debated in the past as it was proposed to be a member of both Agnaridae and Porcellionidae families [3,59].
Our results corroborate the assignment of Mongoloniscus to Agnaridae family as well as the close evolutionary relationship with Porcellionidae family. This relationship is also con rmed by a recently published study undermining the monophyly of Porcellionidae [24].
At family lever our results provide evidence against the monophyly of Agnaridae, Platyarthridae, Trachelipodidae, Porcellionidae and Trichoniscidae families (S3 Fig.). Based on morphological or/and genetic data previous studies came to the same results questioning also the monophyly of Scyphacidae, Philosciidae, Dubioniscidae, Trachelipodidae and Porcellionidae [8,9,17,24]. Aforementioned studies commended on the validity of traditionally used morphological characters for taxonomic purposes and shown the impact of the selected genetic loci on phylogenetic reconstructions.

Conclusions
Present Oniscidea classi cation was determined by morphological characters which probably cannot adequately portray the evolutionary history of the group and its relation with aquatic relatives. This is highlighted by the numerus taxonomic revisions, even at high taxonomic levels, considering new morphological or/and molecular data.
Incongruent patterns were produced using different characters or combination of genetic loci. Herein, keeping separately the genetic data from different organelles, that are not under the same evolutionary pressure, we reached at different phylogenetic patterns even regarding the monophyly of Oniscidea.
Previously published studies indicated the important in uence of taxonomic sampling, outgroup selection, targeted loci and selected models of nucleotide evolution on constructed phylogenies [9,15,23]. Analyses with deviations on any of these parameters could lead in different tree topologies, even between the basal Oniscidea clades, indicating alternative scenarios about the origin and terrestrialization history of the group. Even exploiting similar datasets depending on aforementioned factors we might conclude having different phylogenies [15,23].
Although recently there are growing evidence against the monophyly of Oniscidea [9,15] crucial aspects on taxon's terrestriallity history still remain to be explored. From the phylogenetic point of view, a more inclusive mitogenome dataset could give new insights related to the evolution of both the taxon and the Isopods special mitochondrial structure and organization. A combination of better representation of taxa and nuclear loci should also be applied for more robust and reliable results. In this direction next generation sequencing techniques (NGS) targeting the whole genome, in combination with bioinformatic progress on the development of new advanced tools are anticipated to take advantage of high throughput data and shed light in various facets of terrestrial isopods evolutionary history. Availability of data and materials Genetic data used in the present study are deposited at GenBank and publicly accessible through the provided accession numbers, and other relevant data are presented within paper and Suppl. Materials.

Competing interests
The authors declare that they have no competing interests.