Gene characteristics of the complete mitochondrial genomes of Paratoxodera polyacantha and Toxodera hauseri (Mantodea: Toxoderidae)

Sampling collection and DNA extraction

Two samples P. polyacantha and T. hauseri were collected from Borneo island in 2015, identified by JY Zhang and stored in 100% ethanol at −40 °C. Total DNA was extracted from muscle of one leg using the QIAGEN DNeasy Blood and Tissue Kit (QIAGEN, Germany).

PCR amplification and sequencing

Two mantis mitogenomes were amplified with six pairs of mantis-specific universal primer sets F2, F3, F7, F9, F10 and F11 as described in Zhang et al. (2018a) and specific primers were designed based on the sequenced PCR information from universal primers using Primer Premier 5.0 (Table 1). We used both normal PCR (product length < 3,000 bp) and Long-PCR (product length >3,000 bp) methods with Takara Taq and Takara LATaq DNA polymerase, respectively (Takara, Dalian, China) in a 50 µL reaction volume. The reaction systems and cycling conditions for normal PCR and Long-PCR were as in Zhang et al. (2018a). All PCR products were sequenced in both directions using the primer-walking method and ABI3730XL by Sangon Biotech Company (Shanghai, China).

Mitogenome annotation and sequence analyses

Contiguous sequence fragments were assembled using DNASTAR Package v.6.0 (Burland, 2000). The tRNA genes and their potential cloverleaf structures were identified by MITOS (http://mitos.bioinf.uni-leipzig.de/index.py) (Bernt et al., 2013) using the invertebrate mitogenome genetic code. Two rRNA genes (12S and 16S rRNA) were determined by comparison with homologous sequences of mtDNA from other mantis species using Clustal X (Thompson et al., 1997). Following identification of tRNAs and rRNAs, 13 protein-coding genes were translated with the invertebrate mitogenome genetic code to find the open reading frames between tRNAs (Cameron, 2014b). We used CG View server V 1.0 (Grant & Stothard, 2008) to draw the mitochondrial genome map with GC content and GC skew of P. polyacantha and T. hauseri. The A+T content, codon usage and relative synonymous codon usage (RSCU) of protein-coding genes were calculated by Mega 7.0 (Kumar, Stecher & Tamura, 2016). Composition skewness was calculated according to the following formulas: AT-skew = (A − T)/(A + T); GC-skew = (G − C)/(G + C) (Perna & Kocher, 1995).

Phylogenetic analyses

In order to discuss the phylogenetic relationships of Toxoderidae, 43 previously sequenced mantis mitogenomes (Cameron, Barker & Whiting, 2006; Wang et al., 2016; Ye et al., 2016; Tian et al., 2017; Zhang & Ye, 2017; Zhang et al., 2018a; Zhang et al., 2018b) were used in the phylogenetic analyses. The outgroup taxa were two cockroaches, Cryptocercus kyebangensis and Eupolyphaga sinensis (Zhang et al., 2010) and two termites, Termes hospes (Dietrich & Brune, 2016) and Macrotermes barneyi (Wei et al., 2012). Accession numbers of all mitogenomes are listed in the phylogenetic trees. According to the phylogenetic analyses of Zhang et al. (2018a), we used the nucleotide sequences of the 13 protein-coding genes as the dataset to construct the BI and ML phylogenetic trees. Each of 13 protein-coding genes were aligned using Clustal W in the program Mega 7.0 (Kumar, Stecher & Tamura, 2016) and conserved regions were identified by the program Gblock 0.91b (Castresana, 2000). The resulting alignments were concatenated with Geneious 8.1.6 (Kearse et al., 2012). We used the program PartitionFinder 1.1.1 (Lanfear et al., 2012) to infer the optimal partitioning strategy and choose the best model according to the Bayesian Information Criterion (BIC). The data blocks were defined by each of three codon positions for the thirteen protein-coding genes and a total of 11 partitions were found. ML analysis was implemented in RAxML 8.2.0 with a GTRGAMMA model and branch support for each node was evaluated with 1,000 replicates (Stamatakis, 2014). BI analysis was implemented in MrBayes 3.2 with a GTR + I + G model, each of four chains (three hot and one cold), with run length of 10 million generations and sampling every 1,000 generations (Ronquist et al., 2012). Convergence was assessed with Tracer 1.5 (Rambaut & Drummond, 2007) and trees from the first 25% of the samples were removed as burn-in.

Mitogenome organization and composition

We annotated and deposited the complete mitogenomes of P. polyacantha (MG049920) and T. hauseri (KX434837) in the GenBank database. There two new mitogenomes were double circular DNA molecules with lengths of 15,999 bp and 15,616 bp, respectively (Figs. 1A–1B). These were longer and shorter, respectively, than the mitogenome of S. porioni (15,846 bp), a previously sequenced member of Toxoderidae that we used as a comparison species. The size variation of the three mitochondrial genomes was mainly caused by different intergenic nucleotides (IGNs) and the presence of additional copies of tRNAs in P. polyacantha and S. porioni. Of all 43 sequenced mantis mitogenomes (Cameron, Barker & Whiting, 2006; Wang et al., 2016; Ye et al., 2016; Tian et al., 2017; Zhang & Ye, 2017; Zhang et al., 2018a; Zhang et al., 2018b), the length of the mitogenome of Hierodula patellifera (16,999 bp) was the longest whereas that of Tenodera sinensis (15,531 bp) was the shortest. The mitogenome lengths of seven Paramantini species were long (>16,000 bp) because of a large non-coding region (400–1,500 bp) between trnM and ND2 apart from the typical A+T-rich region. The mitogenomes of Anaxarcha zhengi, Deroplatys desiccate, Mantidae sp., Parablepharis kuhlii asiatica, Phyllothelys sp1., Theopompa sp.-YN and Theopompa sp.-YN were also longer than 16,000 bp because of a long typical A+T-rich region (>1,100 bp) (Ye et al., 2016). The mitogenome of T. hauseri contained the typical 37 genes (13 PCGs, 22 tRNAs and 2 rRNAs) and an A+T-rich region (Table S1) whereas the mitogenome of P. polyacantha had an extra four tRNAs (2 trnA and 2 trnR) (Table S2); by comparison, S. porioni had an extra two trnK. The T. hauseri mitogenome contained the shortest IGNs, a total of 210 bp, compared to 305 bp for P. polyacantha and 313 bp for S. porioni. The nucleotide composition of the P. polyacantha and T. hauseri mitogenomes had a high A+T bias of 74.81% and 73.49% and both showed positive AT-skew and negative GC-skew, which was also similar to S. porioni (Table 2).

Protein-coding genes and codon usages

All 13 protein-coding genes (PCGs) were identified in the mitogenomes of P. polyacantha and T. hauseri. Nine PCGs (ND2, COX1, COX2, ATP8, ATP6, COX3, ND3, ND6 and CYTB) were coded on the majority strand (J-strand) and the remaining four (ND5, ND4, ND4L, and ND1) were coded on the minority strand (N-strand). The length, codon usages and A+T content of PCGs in the P. polyacantha and T. hauseri mitogenomes were nearly identical to S. porioni. Among three mitogenomes, 12 PCGs used ATN (N represents A, T, C, G) as initiation codons with the exception of COX1 which was initiated with TTG. TTG is an accepted conventional initiation codon for many insect mitogenomes including among mantises (Ye et al., 2016; Zhang & Ye, 2017; Zhang et al., 2018a) and cockroaches (Jeon & Park, 2015; Cheng et al., 2016). TAA was commonly used as for the termination codons although the incomplete termination codon T was found in COX3 and ND5 in all three mitogenomes. An incomplete termination codon has also been found in all other sequenced mantis species (Cameron, Barker & Whiting, 2006; Wang et al., 2016; Ye et al., 2016; Tian et al., 2017; Zhang & Ye, 2017; Zhang et al., 2018a; Zhang et al., 2018b). It has been demonstrated that incomplete termination codons can act as functional termination codons in polycistronic transcription cleavage and polyadenylation processes (Ojala, Montoya & Attardi, 1981; Du et al., 2016). In the P. polyacantha mitogenome, COX2 used TAG as the termination codon. Although TAG is the canonical termination codon in insect mitogenomes, it is not used frequently perhaps due to the high percentage of AT nucleotide use by the protein-coding genes (Liu et al., 2016). In the 43 published mantis mitogenomes, only COX1 of Theopompa sp.-YN (Ye et al., 2016) and Leptomantella albella (Wang et al., 2016), COX2 of Theopropus elegans, ATP8 of Sibylla pretiosa, ATP6 of Phyllothelys spp. and Creobroter jiangxiensis, ND4 of Schizocephala bicornis and CYTB of Creobroter urbanus (Zhang et al., 2018a) as well as ND3 of Tamolanica tamolana (Cameron, Barker & Whiting, 2006) used TAG as the termination codon.

The average AT contents of the 13 PCGs in P. polyacantha and T. hauseri were 74.43% and 73.09%, both slightly higher than S. porioni (73%). The PCGs encoded by the majority strand displayed T-skews (the content of T > A) and G-skews (G > C) whereas the minority strand displayed T-skews and C-skews (C > G). We calculated the relative synonymous codon usage (RSCU) of the three Toxoderidae species mitogenomes (Figs. 2A–2C; Table 3) and the result showed that NNU and NNA were higher than 1.0 with the exception of Leu (CUR) and Ser (AGU) in P. polyacantha, T. hauseri and S. porioni and Arg (CGU) only in S. porioni. The most frequent amino acids in the coding sequences of P. polyacantha, T. hauseri and S. porioni mitochondrial proteins were Leu (UUR), Ile and Phe (>300) (Fig. 3). These three amino acids were also frequently used in 43 other mantis mitogenomes (Cameron, Barker & Whiting, 2006; Wang et al., 2016; Ye et al., 2016; Tian et al., 2017; Zhang & Ye, 2017; Zhang et al., 2018a; Zhang et al., 2018b) and Lepidoptera mitogenomes (Liu et al., 2016; Xin et al., 2018; Wu et al., 2012). The least used amino acid in the three mitogenomes was Cys (<50).

Ribosomal RNAs and transfer RNAs

The mitogenomes of P. polyacantha and T. hauseri each had one 16S rRNA and one 12S rRNA gene. The 16S rRNA gene was located between trnL (UUR) and trnV and the 12S rRNA gene was located between trnV and the A+T-rich region as also occurs in other mantises. The size of the 16S rRNA was 1,387 bp in P. polyacantha and 1,322 bp in T. hauseri, both a little longer than in S. porioni (1,317 bp). The size of the 12S rRNA was 787 bp in T. hauseri similar to S. porioni (788 bp) whereas the P. polyacantha (842 bp) was much longer. In the P. polyacantha mitogenome, the A+T content of the rRNA genes was the highest (77.39%) whereas the A+T content of rRNAs in the T. hauseri mitogenome was approximately 76%, slightly lower than the A+T-rich region. In P. polyacantha and T. hauseri, the AT-skew was slightly negative whereas the GC-skew was strongly positive indicating that the contents of T and G were higher than those of A and C, respectively.

Unlike the typical set of 22 tRNA genes in metazoan mitogenomes, there were 26 tRNA genes including an extra two copies of trnA-trnR predicted in the P. polyacantha whereas T. hauseri had a typical set of 22 tRNA genes. The secondary clover-leaf structures of tRNA genes identified in the mitogenome of P. polyacantha and T. hauseri are shown in Figs. 4A–4Z and 5A–5V. The lengths of these tRNA genes varied from 63 bp to 72 bp. All the predicted tRNAs displayed the typical clover-leaf secondary structure, except for trnS (AGN), where the DHU arm appears to be replaced by unpaired nucleotides, a feature typical of other animal mitochondria (Wolstenholme, 1992). All the mismatched base pairs found were U-G pairs, and there were 24 and 25 mismatched base pairs in the P. polyacantha and T. hauseri sequences, respectively. In addition, unmatched U–U base pairs were observed in both two mitogenomes.

The mitogenome of P. polyacantha had an extra 2 copies of trnA-trnR and formed the gene cluster ND3-trnA-trnR-trnA-trnR-trnA-trnR-trnN-trnS-trnE-trnF-ND5. Three trnA genes were identical whereas the first trnR was a little different from the other two because it had an extra 2 bp nucleotide “TA” on the TψC arm (Figs. 4J–4O). Duplication of tRNA is a common phenomenon in mantis mitogenomes. Ye et al. (2016) found that 4 mantises (Creobroter gemmata, Mantis religiosa, Statilia sp., Theopompa sp.-HN) had 2–8 copies of trnR. In addition, Statilia sp. had 5 trnW2 copies as well as 6 trnR forming the gene cluster ND3-trnA-trnR-trnR-trnW2-trnR-trnR-trnW2-trnW2-trnR-trnW2-trnN-trnS-trnE-trnF-ND5. Zhang et al. (2018a) found that the mitogenomes of S. porioni and Schizocephala bicornis had 3 identical trnK and 5 identical trnI, respectively. Evidence for tRNA duplication has also been found in other insect orders such as Ephemeroptera (Zhang et al., 2008), Hymenoptera (Cameron et al., 2008) and Lepidoptera (Hu et al., 2010). However, the occurrence of three copies of trnA-trnR in P. polyacantha is the first such report in Insecta although in the mitogenome of Brontostoma colossus (Kocher et al., 2014), trnA-trnR may have been duplicated at least once. This is because the gene cluster ND3-trnR-trnA-trnR-ND5 was found in B. colossus and there were 40 bp intergenic nucleotides between ND3 and trnR, 20 bp of which showed high similarity to trnA (100%).

A+T-rich region and intergenic regions

The A+T-rich region of the insect mitogenome is equivalent to the control region of vertebrate mitogenomes and harbors the origin sites for transcription and replication (Andrews et al., 1999; Yukuhiro et al., 2002; Cameron, 2014a; Du et al., 2016). The A+T-rich regions of P. polyacantha and T. hauseri mitogenomes were located between 12S rRNA and trnI with lengths of 709 bp and 687 bp, respectively (Figs. 1A–1B), which is similar to S. porioni (685 bp). The length of the A+T-rich region is variable, generally ranging from 639 bp in Mantis religiosa to 1,775 bp in Theopompa sp.-HN (Ye et al., 2016). In the mitogenomes of P. polyacantha and T. hauseri, the contents of A+T were 77.39% and 76.34%, respectively, similar to S. porioni (76.39%). The A +T-rich regions of P. polyacantha and S. porioni genomes both showed positive AT-skew values whereas T. hauseri showed a negative skew. All A+T-rich regions of the three Toxoderidae species displayed negative GC-skew values. Unlike the species Anaxarcha zhengi, Hierodula formosana, Rhombodera valida, Tamolanica tamolana, Theopompa sp.-YN and Theopompa sp.-HN that showed different copies of tandem repetitive sequences (Cameron et al., 2008; Ye et al., 2016; Zhang & Ye, 2017), we failed to find any tandem repetitive sequences in P. polyacantha and T. hauseri using Tandem Repeat Finder V 4.07 (http://tandem.bu.edu/trf/trf.html) (Benson, 1999).

Intergenic regions

The mitogenomes of most insect species seem to be economical (Boore, 1999) although large intergenic regions exist in some species. For example, large intergenic regions located between trnM and ND2 were observed in seven Paramantini species with variable lengths ranging from 296 bp in Tamolanica tamolana to 1,541 bp in Hierodula patellifera (Cameron, Barker & Whiting, 2006; Tian et al., 2017; Zhang & Ye, 2017; Zhang et al., 2018a). In the mitogenome of P. polyacantha, there was a total of 305 bp of intergenic space between genes, of which there were 10 locations of intergenic lengths smaller than 8 bp, four locations of intergenic lengths between 10 bp and 20 bp and five locations of intergenic lengths longer than 20 bp (Table S1). The longest intergenic region was located between COX2 and trnK (76 bp), 28 bp of which showed high similarity to trnK (100%) whereas the other 48 bp showed high similarity to trnD (100%). This gene arrangement can be explained by the tandem duplication/random loss mode (TDRL) (Fig. 6A). Firstly, the region of trnK-trnD was tandem duplicated once. Secondly, the random deletion of a portion of one of the trnK-trnD pairs occurred to form a 76 bp partial “trnK-trnD” residue. The mitogenome of T. hauseri contained a total of 210 bp of intergenic space spread over 21 regions with sizes ranging from 1 to 67 bp (Table S2). The longest intergenic region was located between COX2 and trnK (67 bp) and can also be explained by TDRL mode (Fig. 6B). Accordingly, 22 bp were similar to trnK (100%) and the remaining 45 bp were similar to trnD (100%). Thus, we can infer that trnK-trnD was duplicated at least once and then randomly deleted. This feature was also found in the mitogenome of S. porioni, suggesting that it could be synapomorphic for Toxoderidae, because this characteristic has only been found in Toxoderidae whereas the other families of Mantodea have only 0–2 nucleotides located between COX2 and trnK. In addition to the intergenic region between COX2 and trnK, the mitogenome of S. porioni had another two trnK genes and 2 trnD residues forming the arrangement COII-trnK *-trnD *-trnK-trnD *-trnK-trnD *-trnK-trnD-ATP8 (trnK * and trnD * represent tRNA residues) (Zhang et al., 2018a).

Phylogenetic analyses

The phylogenetic relationships inferred from BI analysis (Fig. 7A) and ML analysis (Fig. 7B) had somewhat different topologies, the BI topology being almost identical to that of Zhang et al. (2018a). In BI analysis, the main topology was as follows: (Amantis nawai (Mantidae) + Mantidae sp.) + ((Liturgusidae + (Leptomantella albella + Iridopterygidae)) + ((Haania sp. + (Schizocephala bicornis + Toxoderidae)) + remaining mantises)). However, in ML analysis, the main topology was (Haania sp. + (Schizocephala bicornis + Toxoderidae)) + (((Amantis nawai + Mantidae sp.) + (Liturgusidae + (Leptomantella albella + Iridopterygidae))) + remaining mantises). The difference was mainly caused by the unstable position of three clades: (Amantis nawai + Mantidae sp.), (Liturgusidae + (Leptomantella albella + Iridopterygidae)) and (Haania sp. + (Schizocephala bicornis + Toxoderidae)). We included three Toxoderidae taxa (P. polyacantha, S. porioni and T. hauseri) and they formed a monophyletic clade which was also the sister clade to Schizocephala bicornis (Mantidae). Toxoderidae as the sister clade to Schizocephala bicornis was also found in Zhang et al. (2018a) but was reported to be the sister clade to Oxyothespinae (Mantidae) in Svenson & Whiting (2009). The species of Oxyothespinae (Mantidae) was not including in this study. To address this discrepancy, future studies should increase the number of species used, especially with samples of Oxyothespinae, Tarachodinae and Amelinae.