Digging deeper: new gene order rearrangements and distinct patterns of codons usage in mitochondrial genomes among shrimps from the Axiidea, Gebiidea and Caridea (Crustacea: Decapoda)

Sample collection

Two species belonging to the infraorder Axiidea (C. ceramica, T. australiensis) and three from Caridea (Macrobrachium bullatum, A. lobidens, Caridina cf. nilotica) were collected from different locations in Australia (Table 1). C. ceramica is represented by two individuals: one from a vouchered specimen (Museum Victoria J40715; GenBank accession number KU350630.1); and the other a previously sequenced sample published incorrectly as T. australiensis (KM501040.1) (Gan et al., 2016). The mitogenome sequence for this latter sample is now registered as C. ceramica under the accession number KU726823.1 and has a genetic similarity of 99.8% to the COI region of a vouchered C. ceramica specimen (Museum Victoria J70519) collected from the same general locality.

To maintain continuity with the original NCBI accession number and species name, the original mitogenome sequence lodged on NCBI for T. australiensis was updated with the newly sequenced mitogenome from a vouchered T. australiensis specimen (Museum Victoria J40711; GenBank accession number: KM501040.2). The accession numbers for each species and associated collecting and identification-related information are detailed in Table 1 and Data S3, including voucher numbers for specimens lodged in Museum Victoria, Melbourne (NMV) and the Museum and Art Gallery of the Northern Territory, Darwin (MAGNT).

Next-generation sequencing and mitogenome assembly

Purification of ethanol-preserved tissue and partial whole genome sequencing (2 × 75 bp for T. australiensis and 2 × 250 bp for others) was performed on the Illumina MiSeq platform as previously described (Gan et al., 2016), after which each mitogenome was assembled with IDBA_UD v.1.1.1 (Peng et al., 2012) and annotated using MITOS (Bernt et al., 2013). Circular mitogenome maps were drawn with BRIG v.0.9.5 (Alikhan et al., 2011). Summary statistics including gene boundaries and length, strand, nucleotide composition, intergenic nucleotides, and number of genes were compiled with MitoPhAST v.1.0 (Tan et al., 2015). Alignment of whole mitogenome sequences and calculation of pair-wise nucleotide identities were performed with SDT v.1.2 (Muhire, Varsani & Martin, 2014).

Gene order analysis

Along with the five mitogenomes sequenced in this study, sequences from 28 other complete mitogenomes from the three infraorders were obtained from NCBI’s RefSeq database (Table 2) for comparative analyses. Arrangements of genes for each of these 33 mitogenomes were compared with all other existing decapod mitogenomes in RefSeq to identify potential novel gene orders unreported by previous studies. Mitogenomes that exhibit gene orders differing from that of the pancrustacean ground pattern (Boore, Lavrov & Brown, 1998) were re-annotated with MITOS (Bernt et al., 2013) to confirm that differences observed are not due to misannotations. Any observed misannotations (e.g., missing genes, incorrect gene boundaries) were corrected before proceeding to further comparative and phylogenetic analyses.

Phylogenetic analysis

Mitogenomes listed in Table 2 were subject to phylogenetic analysis to establish the evolutionary relationships of species within each of the infraorders to provide a framework for interpreting mitogenome gene rearrangements. MitoPhAST v.1.0 (Tan et al., 2015) was used to extract individual PCG amino acid sequences, and these protein sequences were then separately aligned with MAFFT v.7.222 (Katoh & Standley, 2013), followed by trimming with trimAl v.1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). For nucleotide level analyses, PCG nucleotide sequences were manually extracted and fed to TranslatorX v.1.1 (Abascal, Zardoya & Telford, 2010), which aligns nucleotide sequences guided by amino acid translations and then trimmed with Gblocks v.0.19b (Castresana, 2000). On the other hand, rRNA was aligned with MAFFT v.7.222 (mafft-linsi) (Katoh & Standley, 2013) and trimmed with trimAl v.1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). Finally, mitochondrial PCG and rRNA sequences were concatenated into super-alignments to make up the following datasets:

1. 13 PCG (aa) [3,591 characters]

2. 13 PCG (nt) [9,642 characters]

3. 13 PCG (aa) + 12S rRNA + 16S rRNA [5,694 characters]

4. 13 PCG (nt) + 12S rRNA + 16S rRNA [11,755 characters]

Maximum-likelihood (ML) tree inference with ultrafast bootstrap (UFBoot) branch supports (Minh, Nguyen & von Haeseler, 2013) was performed using IQ-TREE v.1.5.0 (Nguyen et al., 2015), which also implements model selection to find the best-fit partitioning scheme. Super-alignments for all datasets were partitioned based on genes. An additional analysis for Dataset B that further partitions it according to first, second and third codon positions was also performed. For Bayesian inference, the same super-alignments generated from all datasets were analysed using Exabayes v.1.4.2 (Aberer, Kobert & Stamatakis, 2014). For each analysis, four independent runs were carried out concurrently for five million iterations each with 25% of initial samples discarded as burn-in. Convergence of chains was checked by ensuring the average standard deviation of split frequencies (asdsf) is below 0.5%, considered to be good convergence according to the Exabayes user guide. Alignments, partitions and best-fit partitioning schemes for all datasets are available as Data S1.

Codon usage analysis

Codon usage (in counts) was calculated using EMBOSS v.6.5.7 (Rice, Longden & Bleasby, 2000) followed by minor adjustments based on the Invertebrate Mitochondrial Code (genetic code = 5). Comparisons among the three lineages of shrimps were made by applying the chi-square test to the pooled codon usage counts for species from each infraorder. Relative synonymous codon usage (RSCU) values were calculated by taking the ratio of the number of times a codon appears to the expected frequency of the codon if all synonymous codons for a same amino acid are used equally (Sharp & Li, 1987). Patterns of variation among individuals in RSCU values were summarized using multidimensional scaling (MDS) based on Euclidean dissimilarities implemented in XLSTAT v.2015.4.01.20978 (Addinsoft, 2010).

Mitogenome composition

Mitogenomes for specimen J40715 of C. ceramica (16,899 bp, 130× cov), T. australiensis (15,485 bp, 86× cov), M. bullatum (15,774 bp, 27× cov), A. lobidens (15,735 bp, 60× cov) and Caridina cf. nilotica (15,497 bp, 63× cov) were assembled into complete circular sequences, annotated (Fig. 1), and deposited in GenBank with accession numbers listed in Table 1. Four of the mitogenomes contain the typical 13 PCG, two ribosomal RNA genes, 22 transfer RNA genes and one long non-coding region while A. lobidens has an additional trnQ flanked by the ND4L and trnT genes (Fig. 1). Detailed composition of each mitogenome can be found in Table S1 while general information on % AT and lengths for all mitogenomes included in this study are in Table S2. The mitogenomes are AT rich (58.9–69.7%), with A. lobidens having the lowest AT content, matching closely to Alpheus distinguendus (60.2%), the only other species of Alpheus having a published mitogenome sequence. Gene lengths are typical but Callianassa and Trypaea have an elevated proportion of intergenic nucleotides, with some spacers in the order of 200 bp in length. This is significantly larger than for other members of the Axiidea, but similar to spacers reported for the Gebiidea. C. ceramica (KU350630.1) has an unusually long control region of 2,036 bp, whereas for all the other taxa it is less than 1,000 bp, including the closely-related T. australiensis (587 bp). However, the control region for the conspecific C. ceramica (KU362925.1) is very similar (1,978 bp) and these two specimens are also very similar in terms of % AT (69.7 and 70.2%). A matrix of pair-wise identities of whole mitogenome sequence alignments can be found in Fig. S1.

Manual inspection and MITOS annotation identified multiple erroneously annotated crustacean mitogenomes in GenBank RefSeq database

The re-annotated species with gene orders divergent to that of the pancrustacean ground patterns obtained from GenBank’s RefSeq database identified several that require revision. One was found to have a missing protein-coding gene and an extra trnS, and others inverted rRNA and tRNA coordinates and other annotation anomalies as detailed in Table 3. For our study, entries were edited based on MITOS annotations and the revised GenBank files for these entries are included as Data S2.

Whole mitogenomes are consistent with monophyly of infraorders Axiidea, Gebiidea, and Caridea

A total of 33 mitogenome sequences from the three groups of interest were utilized (Axiidea: eight, Gebiidea: five, Caridea: 20) to elucidate phylogenetic relationships, with an additional 12 Dendrobranchiata mitogenomes as an outgroup (Table 2). Trees constructed from every dataset and analysis are available in Fig. S2. All trees place M. bullatum, A. lobidens, and Caridina cf. nilotica as sister taxa to other species from their respective genera with relatively high nodal support (UFBoot ≥ 93%, PP 1.00). These trees also consistently place C. ceramica as sister to T. australiensis (UFBoot 100%, PP 1.00) within Axiidea.

The Bayesian-inferred phylogenetic tree in Fig. 2A, constructed based on amino acid sequences of 13 PCGs and nucleotide sequences of two rRNAs (Dataset C), shows a tree topology that is shared by most trees inferred in this study. Most nodes received maximal support from each analysis and dataset. The greatest levels of uncertainty in terms of phylogenetic placement are mostly relating to the relationships among closely related taxa, such as in the Palaemon (6/10 trees), Rimicaris–Opaepele–Alvinocaris–Nautilocaris (2/10 trees), Upogebia–Austinogebia (1/10 trees), Corallianassa–Paraglypturus (1/10 trees), and Macrobrachium (1/10 trees) clades. The only deeper clade with low support is the placement of Atyidae as the sister clade of Alvinocarididae within Caridea (UFBoot ≥ 90, PP 1.00). Within Axiidea, maximal support is observed for almost all nodes. The most basal split in this infraorder separates Strahlaxiidae (Neaxius acanthus) from Callianassidae. Within Callianassidae, the three major lineages correspond to accepted subfamilies, two represented by one species each and Callianassinae by four species. In as far as it goes, the phylogenetic placements of mitogenome sequences in Axiidea are congruent with the current classification at the family, subfamily, and generic levels (Felder & Robles, 2009) (Fig. 2A). These analyses of just four upogebiid species indicate that Upogebia may be paraphyletic with respect to Austinogebia. Since these two genera are nominally represented by over 120 and eight species respectively, any comment on their status is premature at this stage. The degree of divergence between the species of Rimicaris and Opaepele is small relative to the degree of divergence between congeneric species within Macrobrachium, Alpheus, Caridina, and Alvinocaris.

Deviation from the pancrustacean ground pattern is prevalent in the currently sequenced members of the Axiidea and Gebiidea

Most mitogenomes from caridean species have the pancrustacean ground pattern (pattern Gr, Fig. 2B). Those that differ show only minor rearrangements involving the short tRNA genes (patterns Pa, Ap1, and Ap2). In contrast, the mitogenomes of species within Axiidea and Gebiidea exhibit relatively substantial differences in gene order (patterns Up, Ax1, and Ax2) entailing rearrangements of PCG, rRNAs, and a number of tRNAs. An example is pattern Ax2 shown in Fig. 2B, which includes the inversion of the ND1, lrRNA, srRNA, and trnI genes as a block, in addition to the inversion and translocation of trnD from between trnQ and trnM to a position between trnS2 and the putative control region, as well as new placements for trnL and trnV also evident in pattern Ax1.

All gene order novelties relative to the pancrustacean ground pattern are consistent with the relationships depicted by the molecular phylogeny, and in several cases define taxonomic groups at different levels. In Caridea, pattern Pa is common to the three Palaemon species and pattern Ap1 defines the two Alpheus species. Similarly, within Gebiidea, while Thalassina kelanang (family Thalassinidae) has the ground pancrustacean pattern, the other four species are all members of the family Upogebiidae and are united by novel rearrangements involving several tRNA translocations (pattern Up, Fig. 2B). Pattern Ax1, involving the rearrangements of COIII and several tRNAs, is shared among species of Axiidea. The elements of pattern Ax2 that differ from pattern Ax1 support the node joining Callianassa and Trypaea.

Evidence of significant codon usage bias in mitochondrial genomes at the infraorder level

Figure 3 shows there is strong A+T bias in codon usage across the 33 shrimp mitogenomes. RSCU frequencies demonstrate distinct preference for codons with A or T in the third codon position compared to other synonymous codons. Counts of codons used and RSCU values can be found in Table S3. Among the 62 available codons, the four most used codons in all three infraorders are TTT (Phe), TTA (Leu), ATT (Ile), and ATA (Met), all made up of solely A and T nucleotides. Even so, this preference for A+T codons is stronger in Axiidea and Gebiidea mitogenomes and less so in Caridea, most obvious for amino acids Asp (D), His (H), Asn (N), and Tyr (Y). Statistical comparisons show that, for each amino acid, there are significant differences among the three infraorders in the proportions of the codons being used (p-values in Fig. 3). A separate comparison for species with pattern Ax1and Ax2 within the infraorder Axiidea and also with those with the ground pattern (Gr) do not reveal any substantial difference in their codon usage bias (Fig. S3). The MDS plot shows that, for the most part, members of each infraorder cluster together and are largely distinct from the samples from the other infraorders (Fig. 4). A sample of the Gebiidea, T. kelanang, is a maverick, being placed well inside the Caridea cluster and remote from the other members of its infraorder. It also has a very low AT content for this group, being more similar to caridean shrimps. In this context, it is noteworthy that the species is placed as the most basal member of the Gebiidea and is also the only member of the infraorder that has the primitive pancrustacean gene order, which it shares with all other of the members of the Caridea with the exception of some species with minor derived rearrangements (Fig. 2B).